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ABSTRACT 

. We present results from high-resolution three-dimensional simulations of turbulent in- 

' terstellar gas that self-consistently follow its coupled thermal, chemical and dynamical 

O . evolution, with a particular focus on the formation and destruction of H2 and CO. We 

quantify the formation timescales for H2 and CO in physical conditions corresponding 
to those found in nearby giant molecular clouds, and show that both species form 
^ ' rapidly, with chemical timescales that are comparable to the dynamical timescale of 

Ci . the gas. 

We also investigate the spatial distributions of H2 and CO, and how they relate 
fvq . to the underlying gas distribution. We show that H2 is a good tracer of the gas 

^ ' distribution, but that the relationship between CO abundance and gas density is more 

complex. The CO abundance is not well-correlated with either the gas number density 
n or the visual extinction Ay: both have a large influence on the CO abundance, but 
the inhomogeneous nature of the density field produced by the turbulence means 
that n and Ay are only poorly correlated. There is a large scatter in Ay, and hence 
CO abundance, for gas with any particular density, and similarly a large scatter in 
' density and CO abundance for gas with any particular visual extinction. This will 

, have important consequences for the interpretation of the CO emission observed from 

real molecular clouds. 

Finally, we also examine the temperature structure of the simulated gas. We show 
that the molecular gas is not isothermal. Most of it has a temperature in the range of 
10-20 K, but there is also a significant fraction of warmer gas, located in low-extinction 
regions where photoelectric heating remains effective. 

Key words: astrochemistry - molecular processes - methods: numerical - ISM: 
clouds - ISM: molecules 



1 INTRODUCTION 

As e ssentially all star formation within the Milky Way occurs within cold, dense clouds of molecular gas l Lada fc Ladal 



2003), understanding how these clouds form, a nd how and why they then form stars, is crucial if w e are to develop any real 
unde rstanding of the process of star formation ( Mac Low fc Klessen 2004; B allesteros-Paredes et alj|2007t iMcKee fc Ostrikei 



2007). Observations of molecular emission lines provide us with a great deal of information on the internal st ructure and 



dynamics of molecular clouds, and may also hold important clues to their past dynamical history ( Ferrierel 2001 ) 

Traditionally, molecular clouds have been viewed as quasi-static objects that form s tars slowly over a long lifetime 
I Zuckerman fc Evans! 1974 : Shu. Adams fc Lizanol Il987l : iKrumholz. Matzner fc McKeelEoO^ ). n this picture, the dynamical 
evolution of a cloud and the chemical evolution of the gas within it are only loosely coupled, and can be modelled separately. 
In the past few years, however, this traditional picture has begun to give way to a new pictur e of clouds as inherently 
dynamical entities, whose formation and evolution ar e dominated by the effec t s of turbulence jfjallesteros-Paredes et al 



19991 : iKlessen. Heitsch fc Mac Lowlboool : iKlessenl [ioOll : lElmegreen fc Scald liool IScalo fc Elmegreenll2004l ). The dynamical 
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evolution of the cloud is rapid, with a timescale comparable to those of key chemical processes such as the conversion of 
atomic to molecular hydrogen (H2) or the freeze-out of molecules onto the surfaces of interstellar dust grains in dark cores. 
If this picture is correct, then the dynamics and chemistry of the gas are strongly coupled, with one directly influencing the 
evolution of the other, and to model them correctly we must model them together. 

An additional impetus towards the development of coupled models of molecular cloud dynamics and chemistry comes 
from the realization that if the internal motions of the clouds are dominated by supersonic turbulence, then it is far more 
difficult than is often appreciated to infer details of their three-dimensional structure from emission or absorption line obser- 
vations. Because we see clouds in projection, we have no direct information about the line-of-sight positions of emission or 
absorption features, only about their radial velocities. In a supersonically turbulent cloud, random velocity variations along 
the line of sight can create coherent features in position-position-velocity (PPV) space that actually correspond to multiple, 
physically-s eparated regions in real space, or ? conv ersely, ca n break up a real physical feature into multiple distinct features in 
PPV space ( Ballesteros-Paredes fc Mac Low 2002 ; see also Klessen. Heitsch fc Mac Low 200ol . Heitsch. Mac Low, fc Klessen 



2001. Federrath et al 



2009). To properly interpret these observations, we need to be able to compare them with simulated 



observations produced from numerical cloud models that capture the full dynamical and chemical complexity of the clouds 
and the dense cores within them. 

Previous attempts to model cloud formation have not properly addressed the coupling between the cloud dynamics and 
the cloud chemistry. Most attempts have fo cussed either on a detailed treatment of cloud chemistry wi t hin the framework 
of a h ighly simplified dynamical model (e.g. Hennebelle fc Perault 19991 . 200ol ; Kovama fc Inutsuka 2000 , 2002 ; Bergin et al 



20041) or on a detailed treatment of the dynamics, but without any treatment of the chemistry (e.g . Balsar a et al 



2004; 



Kritsuk fc Normanlliool ; Isivz et all [20051 ; Irlennebelle fc Audit! I2OO7I ; jrlennebelle et al.ll2008l ; iBaneriee et al]|2009h . Recently. 



some studies have begun to treat the formation of molecular hydrogen (H2) in high-resolution, three-di m ensional simulation s 



of cloud formation JPobbs. Bonnell. fc Pringlelbood : bobbs fc Bonnellll2007l : Idover fc Mac Lowll2007al rbl: bobbs et al~ll2008l ) 



This is a useful step forward, but although H2 makes up most of the mass of a giant molecular cloud (GMC), it is very 
difficult to observe, owing to the weakness of its rotational lines and their large energy separations. Much of what we know 
about molecular clouds comes instead from observations of carbon monoxide (CO), but until now there has been no attempt 
to model the far more complex chemistry of CO in a high-resolution three-dimensional simulation. 

In this paper, we present a lightweight treatment of gas-phase chemistry and radiative cooling that we have developed to 
tackle this problem. Our treatment is outlined in Section [2] and includes a simplified model for CO formation and destruction 
that tracks the abundances of thirty-two distinct chemic al species f£|2.1l>. an approx imate treatment of H2 self-shielding and 
dust extinction that uses the "six-ray" approximation of Glover fc Mac Low (2007a, b) (i ]2.2|l and a detailed cooling function 
(SJ22J. We show the results of some tests of the model in Section[3] and discuss the initial conditions used for our simulations in 
Section[4] In Sections[5H3 we present a few preliminary results of these simulations, concerning the time evolution of the mean 
chemical abundances (|j5j), the density and temperature probability distribution functions (©, and the spatial distribution 
of molecular hydrogen and CO (iZJ). We conclude with a brief summary in Section [8] 



2 METHOD 



2.1 Basic framework 

We so lve the equations of fluid flow for a magnetised interstellar gas using a modified version of th e ZEUS-MP hydrodynamica l 
code ( Norman 2000;; Haves et al. 2006). An earlier version of this modified code was presented in Glover fc Mac Low (2007a). 
In the present version, we have updated and extended the cooling function, as described in £|2,3l below. We have also improved 
and significantly extended our treatment of gas phase chemistry. We now track the abundances of thirty-two species. Thirteen 



of these species - H 



H 2 , H 3 , CH 



CHJ, OH 4 



H 2 0+, H s O + , CO+, HOC+, O" CT and O^ 



are assumed to be 



instantaneously in chemical equilibrium. For the remaining nineteen species - e , H + , H, H2, He, He + , C, C + , O, + , OH, 
H 2 0, CO, C 2 , O2, HCO+, CH CH 2 a nd CH+ - we follow the full non-equilibrium evolution. 

As in iGlover fc Mac Low! (|2007al ). we represent the abundance of each of the non-equilibrium species with a tracer 
field that advects as a density. To ensure consistent advection of the chemical species, such that elemental abundances are 
conserved locally, as w ell as globally, we use a modified version of the Consistent Multi-fluid Advection (CMA) algorithm of 
Plewa fc Muller ( 19991 ). described in AppendixfX] The chemical rate equations gover ning the creation and destruction of these 
species are solved in an operator-split fashion, using the implicit integrator DVODE ( Brown. Byrne, fc Hindmarshll 19891 ). The 
gas energy equation is also operator-split: the effects of any radiative heating or cooling of the gas are combined with those 
of the pressure-work term into a rate equation that is coupled to the chemical rate equations and so is solved implicitly 
by DVODE, while the advection of the gas energy density is handled as in the unmodified ZEUS-MP code. For reasons of 
computational efficiency, we use conservation laws for charge and elemental abundance to determine the abundances of e~, 
H, He, C and O, reducing the number of rate equations that DVODE must solve to fifteen (fourteen chemical rate equations 
plus the energy equation). 



CO formation in the turbulent ISM 3 



To model the chemistry of the gas, we use a chemical network consisting of 218 reactions between 32 species. Details of 
the reactions included in the network, along with the rate coefficients adopted in each case, are given in Tables IBH4B3I in 
Appendix B. Note that in our present study, we include no grain surface reactions other than the formation of H2 (reaction 
165). The chemical effects of grain surface recombination and the freeze-out of CO, H2O, etc., in dense gas will be treated in 
future work. 



2.2 Photochemistry 

To model the photochemistry of optically th in gas, we ass ume an incident radiation field correspo nding to the st andard 
interstellar radiation field, as determined by Drainei ( 19781 ). This field has a strength Go = 1.7 in Habin d 1 19681 ) units, 



corresponding to an integrated flux of 2.7 x 10 -3 erg cm -2 s _1 . Photodissociation and photoionization rates appropriate for 
this field strength are listed in Table [B2l 

In optically thick gas, it is necessary to account for absorption by dust, by the Lyman- Werner lines of H2 (important for H2 
self-shielding and CO shielding), and by the ultraviolet absorption lines of CO (important for CO self-shielding). To treat these 
processes with full accuracy in a three-dimensional hydrodynamic simulation is beyond our current capabilities. The problem 
is one of computational cost: in a hydrodynamical simulation with N fluid elements, the cost of resolving a monochromatic 
radiation field with the same spatial resolution, and with a comparable angular resolution, is of order 0(N 5 ^ 3 ), as we require 
0(N 2 ^ S ) rays to fully sample the angular distribution of the radiation field for each of the N fluid elements. To model line 
absorption, the radiation field must also be discretised in frequency-space, increasing the cost to 0{N V x N 5 ^ s ), where N v is 
the number of frequency bins required. In comparison, the cost of modelling the hydrodynamics is only of order O(N). Thus, 
the cost of properly treating the radiation field in an optically thick 128 3 zone hydrodynamical simulation is 128 2 iV„ times 
greater than the cost of solving for the hydrodynamical evolution, which is far out of the reach of our current computational 
resources. 

To overcome thi s problem, we are forced t o approximate. The approach that we h ave chosen to adopt is the "six- 
ray" method used in Glover fc Mac Low ( 2007al lbh . which is based on an original idea by Nelson fc Langer (1997). In this 



approximation, we compute photochemical rates in each zone in our simulation volume by averaging over a small number of 
lines of sight. Specifically, we compute the column densities of H2, CO, and H nuclei in all forms (i.e. iVH.tot — A r H+2A r H 2 +^h+ > 
where Nu, Nh_ 2 and N H + are the column densities of atomic H, molecular H2, and protons, respectively) in both the positive 
and negative directions along each of the three coordinate axes of the simulation. For each zone we therefore know the column 
densities in six different directions (or, alternatively, along six rays). 

The rates of most of the reactions listed in Table IB2I are sensitive only to the amount of dust extinction, and for a 
plane-parallel slab geometry, the rates in optically thick gas are related to those in optically thin gas by the expression 

-Rthick = -Rthin/dust = -Rthin exp(— 7t4 V ), (1) 

where ilthin is the optically thin rate, /dust = exp(-^Ay) is the dust shielding factor, Ay is the visual extinction in magnitudes, 
and where the appropriate value of 7 for each reaction is listed in Table lB2l Using the relationship 

. _ iVH.tot ,„,. 

Av ~ 1.87 x 10 21 cm-2 [I > 



between Ay and iVH.tot , as is appropriate for gas in the diffuse ISM (jPraine fc Bertoldil ll996T ) . we can associate a value of 



Ay with each of our six rays for any given zone in our simulation volume. We can then calculate the total rate in that zone, 

-Rthick.i, as 

1 6 

#thick,i = -iithin y^exp(-7Ayj), (3) 



6 



where Ay t j is the visual extinction along ray j. 

To treat the photodissociation of HbO" 1 " and HaO" 1 " (reactions 188-195), we use a very similar appro ach. However, for these 



react ions the optically thick and optically thin rates are related by the more complicated expression (jSternberg fc Dalgarno 
19951 ). 



ilthick = i? thin exp(-2.554v + Cl.0165.4v), (4) 
for Ay ^ 15, and by Equation [T] with 7 = 2.8 for Ay > 15. Hence, the total rate in zone i is given in this case by 

6 

i?thick,i = ^-Rthin [/ A exp(-2.55Av, J + 0.0165A v ,i) + (1 - /a) eacp(-2.8Av,i)] , (5) 
j=i 

where /a = 1 for Ayj ^ 15, and /a = for Ay,j > 15, and where Ay,j is computed in the same manner as before. 

To treat H2 photodissociation accurately, one must take into account not only absorption by dust, but also the effects 
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of H2 self-shielding. If we assume that the effects of dust absorption and self-shielding can be treated separately, then for a 
plane-par allel geometry the e ffects of dust absorption can be modeled by a dust shielding factor given b y Equation [T] with 
7 = 3 .74 ( Draine fc Bertoldi 19961 ). Similarly, the effects of H2 self-shielding can be modeled by a factor ( Draine fc Bertoldi 
19961 ) 



/shield — 



0.965 



0.035 



-^exp [-8.5 x 10~ 4 (1 



\V21 



(6) 



(l + z/b 5 ) 2 (l + x)V 

where x = N n J(5 x 10 14 cm" 2 ), and 65 = 6/(i0 5 cms ), where b is the Doppler broadening parameter. The fully shielded 
H2 photodissociation rate then follows as 

^H 2 ,thick = /dust/shield-RH 2 ,thin- (7) 

Using our six-ray approximation, we can compute /dust and /shield for each ray, and hence can compute the total rate as 



Rn 2 , thick, i — — J?H 2 ,thin ^ ] /dust, j /shield, j ■ 
3=1 



(8) 



Finally, to treat CO photodissociation, it is necessary to take account of three separate contributions to the shielding: CO 
self-shielding, shielding of CO by the H2 Lyman- Werner lines, and dust shielding. For a plane-parallel geometry, the shielded 
CO photodissociation rate is related to the optically thin rate by 

-RCO, thick = /co/H 2 /dustficO,thin, (9) 

where fco and /h 2 are functions of the CO and H2 column densities, respectively, and have been tabulated by Lee et all 
(1996), and where /dust is given by Equation [T] using the value of 7 from Table [B2l Using our six-ray approximation, we can 
compute values of fco , /h 2 and Ay for each ray, from which the partial contribution from that ray to the total rate follows 
as 



-RcO, thick, i = - -BcO.thin /cO,j /h 2 , j /dust ,j ■ 
3=1 



(10) 



2.3 Thermal model 

We model the radiative and chemical heating and cooling of the gas with a cooling function that contain s contributions from 
18 diffe rent processes, listed in Table [1] Our treatment of most of these processes largely follows that in lGlover fc Mac Low 
(2007a), to which we refer readers desiring further details; the few exceptions are noted below. 



2.3.1 Fine structure cooling 



The extremely simplified model of ISM chemistry presented in iGlover fc Mac Low! (|2007al ) assumed that all of the carbon in 
the gas would be kept in singly ionized form by the interstellar radiation field, and so in that paper there was no need to 
treat cooling from the fine structure lines of neutral carbon. However, the significantly improved treatment of ISM chemistry 
presented in this paper removes this assump tion and so it is necessar y to include them in our thermal model. To do this, 
we largely follow the same prescription as in IGlover fc Jappsenl (|2007l ) . The one excepti on is in the rates for the collisional 
excitation of the fine structure lines o f atomic carbon by col l isions with atomic hydrogen. IGlover fc Jappsenl (|2007l ) use rates 
for this process that were taken from Irlollenbach fc McKeel dl989l ). but in this work we use instead the more accurate rates 
recently calculated by Abrahamsson. KremTX*I3algarnor ( 20071 ) . We also adopt their new rates for the excitation of atomic 
oxygen by atomic hydrogen in place of the rates used in our previous work. 



2.3.2 CO and H 2 cooling 

To treat rotational cooling from CO and H2O, we use the tabulated cooling functions of Neufeld fc Kaufman (1993) and 
Neufeld. Lepp fc Melnickl (|1995l ). They use a large velocity gradient (LVG) approach to compute the cooling rates from CO 
and H2O as a function of temperature, density and effective optical depth. Based on the results of their calculations, they 
present fits to the CO and H2O rotational cooling functions of the form 

^l/2-t-O 



1 



nH 2 J_ 

Lute Lq 



n H2 



_™l/2 



(11) 



where L m is a cooling rate coefficient defined such that the cooling rate per unit volume from species m (where m = CO or 
H2O) is given by A m = L m n m nH 2 , and where Lq is the cooling rate coefficient in the low density limit, Llte is the cooling rate 



CO formation in the turbulent ISM 5 



Process 



Reference(s) 



Cooling: 

C fine structure lines 



C+ fine structure lines 



O fine structure lines 



H2 rovibrational lines 
CO and H2O rovibrational lines 
OH rotational lines 
Gas-grain energy transfer 
Recombination on grains 
Atomic resonance lines 
H collisional ionisation 
H2 collisional dissociation 
Compton cooling 



Atomic data - Silva fc Viegas (2002) 
Collisional rates (H) - Abrahamsson, Krcms 
Collisional rates (H2) — ^clvrode^et^l^ (19SH) 
Collisional rates (e _ ) — ^ohnso^e^il^ (1987) 
Collisional rates (H+) - Roueff fe Le Bourlot (1990) 
Atomic data - Silva fc Viegas (2002) 
Collisional rates (H2) — ^lower^; Launay (1977) 
Collisional rates (H, T < 2000 K) - Hollenbach fc McKee 
Collisional rates (H, T > 2000 K) - Keenan et al. (1986) 
Collisional rates (e~) - Wilson fc Bell (2002) 
Atomic data - Silva fc Viegas (2002) 

Collisional rates (H) - Abrahamsson, Krcms & Dalgarno < 
Collisional rates (H2) 
Collisional rates (e~) 

Collisional rates (H+) - Pequignot (1990. 1996) 
Le Bourlot. Pineau des Forets fc Flower (1999) 



see Glover fc Jappsen (2007) 
Bell. Berrington fc Thomas (1998) 



Neufeld fc Kaufman (1993) ; Neufeld. Lepp fc Melnick (1995) 

Pavlovski et al. (20021 

Hollenbach fc McKee (1989) 

Wolfire et al. (2003) 

Sutherland fc Dopita (1993) 

Abel et al. (1997) 



See Table IBTI 
Cen (1992) 



Heating: 

Photoelectric effect 

H2 photodissociation 

UV pumping of H2 

H2 formation on dust grains 

Cosmic ray ionisation 



Bakes fc Tielcns (1994); Wolfire et al. (2003 



Black & Dalgarno (1977) 

Burton. Hollenbach fc Tielens (1990) 

Hollenbach fc McKee (1989) 



Goldsmith & Langcr (1978) 



Table 1. Processes included in our thermal model. 



per molecule when the rotational level populations are in local thermodynamic equilibrium (LTE), and ni/2 is the H2 number 
density at which L m = 0.5Lo. Lo is purely a function of temperature, but t he other three fit parameters (L ute, ^1/2 and 
a) depend on both the temperature and the effective optical depth of the gas. iNeufeld. Lepp fc Melnick i|l995l ) parameterise 
the latter in terms of an effective column density pe r unit velocity, A^(m) , for ea ch coolant m. For a flow without any special 
symmetry in which the LVG approximation applies, Neufeld fc Kaufman ( 19931 ) give this parameter as 



Mm) = 



(12) 



For CO rotational cooling, INeufeld fc Kaufman |l993l) and INeufeld. Lepp fc Melnickl (119951 ) tabulate values of the fitting 
parameters for temperatures in the range 10 K < T < 2000 K and effective column densities in the range 14.5 < logiV(CO) < 
19.0, where iV(CO) has units of cm -2 per km s -1 . For H2O rotational cooling, they tabulat e values for temperatures 1 K < 
T < 4000 K and optical depth parameters 10.0 < log A^(H 2 0) < 19.0. At low temperatures. INeufeld. Lepp fc Melnickl (jl995T ) 
list values of the cooling rate for both the ortho and the para variants of H2O. To compute the total water cooling rate, we 
assume that the ortho:para ratio is fixed at 3:1. 



To properly treat gas below 10 K, we would not only have to extend the lNeufeld. Lepp fc Melnickl (|1995l ) cooling functions 
to lower temperatures, but would also have to take into account several other physical processes that are not included in our 
current model, such as the freeze-out of CO and wate r onto dust grains, or the fact that the dust temperature in the cloud will 
decrease as the extinction increases I Goldsmith! 2001 ) . To avoid this, in the simulations presented here we have introduced an 
artificial temperature floor at 10 K and switch off radiative cooling in gas colder than this. 

To account for cooling from CO and water at very high temperatures (T > 2000 K for CO, T > 4000 K for H2O), we 
adopt cooling rates that are the same as the rates at the largest tabulated temperature. As we expect CO and H2O to be 
rapidly destroyed by collisional dissociation and chemical reactions with atomic hydrogen at these high temperatures, any 
uncertainty in the cooling rates per CO or H2O molecule in this regime will not have a large effect on the total cooling rate, 
owing to the small molecular abundances. This approxima tion should therefore be re asona ble for most uses. Nevertheless, a 
treatment of CO and H2O cooling along the same lines as INeufeld fc Kaufman! jl993h and INeufeld. Lepp fc Melnickl (|l995l ), 
but which extended to higher temperatures, would clearly be desirable. 
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To treat gas with N below the tabulated range, we simply adopt the same fitting parameters as are given for the smallest 
value of N that is tabulated. As the latter generally corresponds to gas that is already very close to the optically thin limit, 
this assumption should give accurate results. To handle the case where N in the simulation exceeds the largest tabulated 
value, we again use rates corresponding to the largest value that is tabulated. Consequently, we will overestimate the cooling 
rate in very dense, highly shielded gas, particularly when the velocity divergence of this gas is small. However, as this gas is 
unlikely to be well-resolved in our simulations in any case, this simplification is again unlikely to introduce large uncertainties 
into our results. 

The lNeufeld fc Kaufman] jl993h and lNeufeld. Lepp fc Melnickl (|l995T ) treatments assume that only collisions with H2 are 
important in determining the CO or H2O rotational cooling rates, as is appropriate in a fully molecular gas. However, at 
early times or at low gas densities, the H2 abundan ce is small, and collisions with at omic hydrogen or with electrons may also 
become important. We therefore follow lYanl (|1997l ) and iMeiierink fc Spaansl §005) and replace nn 2 in Equation [TT] with an 
effective number density n e g. For CO rotational cooling, n e g is given by 



w c ff,co,rot = nu_. 



n H + 



1.3 x 10" 



(13) 



where <r H = 2.3 x 10 15 cm 2 , cr H; 
rotational cooling, n c g is given by 



Weft, H 2 0, rot = n H 2 + 10 Tin 



where fc Ha = 7.4 x lO'^T 1 



fee 

fen. 



= 3.3 x 10" 16 (T/1000 K)" 1/4 cm" 2 , and v e = 1.03 x 10 4 (T/l K) 1/2 cm s" 1 . For H 2 



(14) 



60.191/T 



2/31 



(15) 



a 3 s 1 , and k B is given by 

k c = dex [-8.020 + 15.749/T 1/6 - 47.137/T 1/3 + 76.648/T 1/2 

These for mulae for n e g are adapted from those in IMeiierink fc Spaansl |2005h , with one exception: the exp ression for k c is 
taken from Faure. Gorfinkiel fc Tennyson ( 2004 ). because the expression given by Meiierink fc Spaansl 1 2005h for k e blows up 

at low temperatures. 

To treat CO and H2O vibrational cooling, we again use the results of iNeufeld fc Kaufman! (jl993h . They present fitting 
functions of the form 

Lm 



1 

Lo Lute 



(16) 



for both CO and H 2 0, and give analytical fits to Lo(T) for both coolants, as well as numerical values for Lute for temperatures 
100K < T < 4000K and effective column densities 13.0 < logiV < 20.0. Our treatment of cooling at temperatures and effective 
column densities that lie outside these bounds is the same as that used for CO and H2O rotational cooling, as described above. 
As before, we account for CO and H2O cooling in gas that is not f ully molecular by replacing tih 2 in Equation [16] with an 
effective number density n e g, taken from IMeiierink fc Spaansl (|2005l ). For CO vibrational cooling, this is given by 

I Lco,c 



™eff,CO,vib = ™H 2 + 50 nH + 



\ico,o 



where 

icO.e 

Lco,o 



1.03 x 10" 



V300 7 



1.14 x 10" exp I 



sxp 

/-<is.0\ / 



-3080 \ 
T ) 



-3080 \ 
T ) 



while for H 2 vibrational cooling we have 



n H2 + 10n H + 



n C rT,H 2 0,vib 

where 

£h 2 o, c = 2.6 x 10 



J H 2 0,o 
^H 2 O,0 



-1/2 



exp 



-2325 \ 
T ) 



->H 2 O,0 



0.64 x 10" exp 



(l^) ex p(- 



-2325 \ 



(17) 

(18) 
(19) 

(20) 

(21) 
(22) 



JU/3 j —r y rp 

and where in all of these expressions, T is the gas temperature in K. 

Finally, we note that when iV is large, isotopic variants of CO and H2O can contribute significantly to the total cooling 
rate, as they will have much smaller effective column densities a nd so will be less affected by self- absorption. To model cooling 
from the isotopic species 13 C ie O, 12 C ls O and H 2 ls O, we follow INeufeld. Lepp fc Melnickl |l995l ) and assume that the cooling 



rate coefficients for these isotopic species are the same as for the standard variants 12 C le O and H2 ie O, and that they have 
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abundance ratios of 1%, 0.2% and 0.2%, respectively, compared to the standard variants. We do not include the effects of 
other isotopic variants (e.g. deuterated water, HDO or D2O) as we expect their abundances to be too small for them to 
contribute significantly. 



2.3.3 OH cooling 

To model cooling from OH, we use a rate taken from lPavlovski et al. 1 20021 ) . based on iHollenbach fc McKeel \ 19791 ). This rate 
assumes that the OH molecu les are not in LTE, which is a reasonable assumption provided that the gas density n < 10 10 cm -3 
I Hollenbach fc McKeel Il979h . 



3 CODE TESTS 

In testing our modified version of ZEUS-MP, our main focus was on verifying the ad ditional physics that we h ave added to 



the c ode, as the unmodified code has already undergone significant testing (see e.g. Stone fc Norman 1992al lbt lHaves et al 
20061 ) 



To verify that our simplified model of CO formation and destruction performs as expected, we have used our chemical 
network and cooling function to model the chemical and thermal evolution of static gas for a range of different number densities 
and obscuring average column densities (or visual extinctions). We then compared the results of these single zone models 
with the results of similar simulations performed us ing a detailed chemical network derived from the UMIST99 compilation 
of reaction rates ( Le Teuff. Millar fc Markwick 2000), and using the same cooling function. By modelling the thermal as well 



as the chemical evolution of the gas, we are able to quantify the effect of errors in the chemical abundances of the major 
coolants on the thermal state of the gas. In all of these tests, we use the same incident radiation field, cosmic ray ionization 
rate, elemental abundance of carbon and oxygen, etc. as in our three-dimensional simulations described below. We ran each 
of our test models for a total time of 3.1 x 10 14 s, or about 10 million years. 

Our initial comparisons showed significant discrepancies between the results obtained using our simplified chemical model 
and the UMIST99-derived model. The abundances of the dominant carbon and oxygen carrying species were generally re- 
produced accurately in both models, but we found large differences in the abundances of some of the trace species. However, 
further investigation showed that these discrepancies were caused primarily by differences in the reaction rate coefficients 
adopted for several key reactions in our model compared to those used in the UMIST model. Specifically, we found that to get 
good agreement between the models, it was necessary to ensure that the same values were used for the rates of the following 
reactions: H + recombination (reaction 12); He + recombination (reaction 17); charge transfer from He + to H (reaction 18); 
the destruction of CH by atomic hydrogen (reaction 35); the formation of CO from C and OH (reaction 46); the formation of 
O2 from O and OH (reaction 47); dissociative charge transfer from He + to CO (reaction 104); H^ dissociative recombination 
(reactions 110-111); CH + dissociative recombination (reaction 112); H20 + dissociative recombination (reactions 120-122); 
H3O 4 " dissociative recombination (reactions 123-126); H2 formation on dust grains (reaction 165), and, lastly, H2 photodis- 
sociation (reaction 168). In a couple of cases, the differences in rate coefficients are a result of the use of a different low 
temperature extrapolation from the same experimental data, but in most cases, the differences are due to our use of recent 
experimental or theoretical values that post-date the construction of the UMIST99 model. We therefore expect the values used 
in our model (and listed in Table IBT) to be the more accurate ones. We also note that we are not the first group to remark on 
the sensitivity of astrochemical reaction netwo rks to the large uncertainties that exist in some reaction rate coefficients (see 



Vasvunin et~al1l2004l : IWakelam et~al1l2005l ) , and that uncertainties of this kind in the input physics are an important, but 



e.g. 

currently unavoidable, limitation on the accuracy of simulations of molecular cloud chemistry. 

Having made these adjustments, so that we are comparing like with like, we find good agreement between the results 
produced using our simplified chemical network and the results produced using the full UMIST network. This is illustrated 
in Figures HHS] In Figure [TJl, we show how the C + , C and CO abundances vary as a function of Ay at the end of our test 
runs for an initial density of no = 100 cm -3 . The results from our simplified model are plotted as dotted lines, while those 
from the UMIST model are plotted as solid lines. Figure QJ> gives a similar comparison of the OH, H2O and O2 abundances, 
while in Figure [TJ; we plot the ratio of the gas temperature in the simplified model to that in the full model. We find very 
good agreement for all of the plotted abundances at all Ay, and only very small differences in the temperature. 

In Figure(2] we give a similar comparison for the case of gas with no = 1000 cm" 3 . Again we find good agreement for most 
species, although in this case the simplified model predicts too large an abundance of neutral carbon at Ay > 6. However, 
this disagreement is possibly a little misleading. As Figure [3] demonstrates, the evolution of the C abundance with time at 
high Ay is very similar in both the simplified and the full models, but the final conversion of the residual C and C + to CO 
occurs slightly later in the simplified model. It should also be noted that we expect our neglect of the effects of freeze-out to 
have a larger impact on the gas-phase chemistry of high Ay gas than the small discrepancy noted here. 



Finally, regarding the cooling function, we note that most of its features are already well-tested, as described in ldover fc Mac Low 
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Figure 1. (a) Abundances of C + , C and CO, plotted as a function of Ay, at the end of our static, single-zone simulations with 
no = 100 cm -3 . The results produced by our simplified chemical model are given as dotted lines, while the results of the UMIST model 
are shown by solid lines, (b) As (a), but for the OH, H2O and O2 abundances, (c) As (a), but showing the ratio of gas temperature 
produced by the simplified model to that produced by the full UMIST model. Note that we plot the ratio rather than the individual 
temperatures because the difference between the models is very small, of the order of 0.05%. 



(2007a|). The main addition that requires testing is our implementation of CO and H2O cooling. We have verified that this is 
i mplemented correctly by ens uring that we can reproduce all of th e tabulated values for the CO and H2O cooling rates given 



Neufeld fc Kaufman (|l993h and lNeufeld. Lepp fc Melnickl (jl995l l. and by ensuring that the rates vary smoothly in between 



111 



the tabulated values. 



4 INITIAL CONDITIONS 



In this paper, we present the results of a small set of simulations of supersonic turbulence designed to address the issue of 
numerical convergence, and to highlight the capabilities of the code. We performed three simulations with numerical resolutions 
of 64 3 , 128 3 and 256 3 zones, respectively, that we will hereafter denote as runs Rl, R2, and R3. All three simulations shared 
the same set of initial conditions: a periodic box of side length L = 20 pc, filled with initially uniform atomic gas with a 
density no = 300 cm -3 and a temperature To = 60 K, permeated by a magnetic field with an initial field strength Bo = 2 fiG 
oriented parallel to the z-axis of the b ox. The initial turbu l ent vel ocity field had an RMS velocity of « rm8 = 5 kms" 1 and was 
constructed in the same fashion as in Idover fc Mac Low! d2007bl ). The turbulence was driven as outlined inlMac Low et al 



1 19981 ) and iMac Lo^l (|l999h with a driving power E — 2.805 x 10 erg s so that the RMS velocity of the turbulence 



remained approximately 5 km s _1 throughout each of the runs. To allow us to make meaningful comparisons between the 
different r esolution runs, we e n sured that the same turbulence driving pattern was used in each case. 

As in Idover fc Mac Low! (|2007alibh . we adopted standard solar abundances of hydrogen and helium, and abundances of 
carbon and oxygen taken from Sembach et al. (2000), i.e. xq = 1.41 x 10 -4 and xo = 3.16 x 10 -4 , where xc and xo are the 
fractional abundances by number of carbon and oxygen relative to hydrogen. 
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Figure 2. (a) As Figure[TJi, but for no = 1000 cm" 3 , (b) As (a), but for the OH, H2O and O2 abundances, (c) As (a), but for the gas 
temperature T. 
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We ran each simulation until a time f cn( j = 1.8 x 10 s ~ 5.7 Myr. This corresponds to almost three turbulent crossing 
times, tcross = Z//(2« rms ) ~ 2 Myr. As we shall see later, this is long enough to allow most of the fluid quantities to reach a 
state of statistical equilibrium. 

It should be noted that our simulations are somewhat inconsistent in that they adopt periodic boundary conditions for 
the gas, but do not do the same for the radiation. This is an unfortunate but necessary compromise. Our simulations do not 
have sufficient dynamical range to follow both the process of cloud assembly and the evolution of gas within the cloud, and so 
we use periodic boundary conditions for the gas in order to be able to focus on a small volume of already assembled material, 
while continuing to use a more physically appropriate boundary condition for the radiation that has it simply penetrating 
inwards from the edges of the box. 



5 TIME EVOLUTION OF MEAN CHEMICAL ABUNDANCES 

We begin our discussion of the results of our simulations by examining the time evolution of the spatially-averaged mass- 
weighted mean abundances of several key chemical species in our three simulations. We define the mass-weighted mean 
abundance of a species m as 

, , E I . J . fc :r n 1 (j,J,fc)pA\/(i,i,fc) 

{x m )M = t- , (16) 

Mtot 

where x m (i, j, k) is the fractional abundance of species m, p(i,j, fc) is the mass density in zone (i, j, k), AV(i, j, k) is the volume 
of zone k), M to t is the total mass of gas present in the simulation, and where we sum over all grid zones. In all but one 
case, we define the fractional abundance x m of species m as the abundance by number of the species relative to the abundance 
of hydrogen nuclei, i.e. 

^ m = — , (24) 
n 

where n m is the number density of species m and n is the number density of hydrogen nuclei. The exceptional case is that 
of H2, for which this convention would give a value of xn 2 = 0.5 for gas in which the hydrogen is fully molecular. This has 
the potential to be confusing for readers who are unfamiliar with this convention, and so for clarity we define the fractional 
abundance of H2 to be 

zh 2 = -, (25) 



so that ih 2 = 1 for fully molecular hydrogen. Note also that this is the same convention as is used in iGlover fc Mac Low 



( 2007al ibh. and so the values of xh 2 and (i H2 )m discussed here can be directly compared to those in our previous papers 



5.1 Molecular hydrogen 



In Figuref?] we show how the mass-weighted mean abundance of H2, {x-h 2 )m, evolves in runs Rl, R2, and R3. At early times, 
we see some dependence on the numerical resolution of the simulation. This is probably a consequence of the fact that denser 
structur es are formed as we in crease the numerical resolution, on account of the better resolution of turbulent compression 
(see e.g. iFederrath et al.l 12003 . Figure 5). As the H2 formation rate depends on the density, the formation rate thus also 
increases. 



Similar results were found previously bv lGlover fc Mac Low! (|2007bl ). At late times, most of the resolution dependence 
disappears. The value of (:th 2 }m in the 64 3 zone run remains slightly smaller than that in the higher resolution runs, but 
there is almost no difference between the results of the 128 3 and 256 3 zone runs, suggesting that in this case, a numerical 
resolution of 128 3 zones is enough to reproduce the final H2 abundance accurately. 

Regard ing the time evolution of th e H2 fraction, we first note that we see the same rapid growth in the H2 abundance as 
we found in Idover fc Mac Low! (|2007bh . Within only 1 Myr, the hydrogen has already become 50% molecular. Nevertheless, 
it is also clear that (s;h 2 )m has yet to settle into a steady state by the end of the simulations at t = 5.7 Myr. Although the H2 
chemistry in the denser gas has largely reached a steady state by this point, there remain some low-density regions in which 
the H2 formation timescale is longer than the time elapsed in the simulation (despite the acceleration of H2 formation caused 
by the turbulence). 



5.2 Carbon chemistry: C , C and CO 

In Figure [5] we examine the time evolution of the mass- weighted mean abundances of C + , C, and CO. The first point to 
note here is the very rapid change in (:e c +)m and {xc)m at the beginning of the simulation. This is a consequence of our 
choice of initial conditions. We began with all of our carbon in the form of C + , which has a short recombination timescale 
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^rec,c+ < 0.1 Myr for our initial temperature of 60 K and initial density of 300 cm -3 . In the interior of the simulation volume, 
the value of x c + in photoionisation equilibrium is typically much smaller than this starting value. Thus, there is a rapid 
conversion of C + to C, resulting in the rapid changes in {x c +)m and (ic)m visible in the Figure. Once we are past this initial 
transient, the subsequent evolution of {x c +)m and (xc)m occurs on a much slower timescale. It is driven by a combination 
of two main factors: the changes that are occuring in the density structure of the gas, in response to the turbulence, and the 
conversion of C + and C into CO. 

The CO abundance evolves rapidly within the first few million years of the simulations. During the first 0.5 Myr, there 
is almost no CO present, but by t = 1 Myr, 12% of the total carbon has been incorporated into CO. This fraction increases 
to 32% by t = 2 Myr, and reaches a steady-state value of around 50% at t > 4 Myr. Atomic carbon dominates for t < 3 Myr, 
while CO dominates at t > 3 Myr. The steady-state valu e of the C/CO ratio is approximately 55%. This is consistent with 
the values computed by IPapadopoulos. Thi fc Vitil ( 2004 ) for gas with a density of order 10 3 cm -3 and a visual extinction 
Ay ~ 1-10 that is illuminated by the stan dard interstellar radiation fi eld, and is compatible with the ratio observed in a 
number of nearby molecular clouds (see e.g. iPlume. Jaffe fc Keend ll994). 

At early times, the C and CO abundances clearly depend on the numerical resolution of the simulation: in higher resolution 
simulations, we find more CO and less C than in lower resolution simulations. However, the difference between the runs lessens 
with time and with increasing resolution. For times t > 4 Myr, there is very little difference between the results of the 128 3 
and 256 3 zone runs. The C + abundance shows very little sensitivity to the numerical resolution throughout the simulation, 
as it is located primarily in large, low-density regions that are well resolved in all of our simulations. 



5.3 Oxygen chemistry: O, OH, H2O and O2 

In Figure [6] we examine the time evolution of the mass-weighted mean abundances of O, OH, H2O, and O2. It is clear from 
the Figure that most of the oxygen remains in atomic form at the end of the simulation, with roughly 23% having been 
incorporated into CO (not plotted here), and less than 1% into other molecules (primarily H2O at early times and O2 at 
late times). Our values for the O and OH abundances appear to be numerically well converged during most of the period 
simulated, with the exception of a short period around t = 2 Myr; note that this corresponds to a turbulent crossing time, 
and hence to the time at which the first major shock-shock interactions are occurring. This is also the time at which most 
of the CO is forming. Our values for the H2O and O2 abundances are less well converged, although there is some indication 
that the results of the 128 3 and 256 3 runs have converged by the end of the simulation. The O and OH abundances appear 
to have reached a steady state by t — t cn d, but there is no indication that the H2O and O2 abundances have done so. 

If we compare the mean mass- weighted abundances of H2O and O2 that we find in these simulations with the values 
(or upper limits) measured in local star-forming regions, then it becomes immediately apparent that there is a significant 
discrepancy. In our simulations, we find that at times t > 2 Myr (corresponding to roughly a single turbulent crossing time), 
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(sh 2 o)m ~ 3-5 x 10 -7 and (a;o 2 )M ~ 1-10 x 10 -7 . However, observations of a number of local star-forming regions with 
the Submillimeter Wave Astronomy Satellite (SWAS) and the Odin satellite find H2O and O2 abundances that are more 
than a factor of ten smaller (see e.g. Bergin et al. 2000l : Goldsmith et al. 2000l ; Pagani et al. 20031 ; Larsson et al. 2007 ). This 
discrepancy is probably due to the neglect of freeze-out processes in our current study. Static gas-phase chemical models of 
mole cular clouds overproduce H2O and O2 in a similar fashion to our dynamical models ( Bergin et al. 2000l : Goldsmith et al 



2000) and the inclusion of grain-surface processes in these models is widely seen as the most promising way to restore agreement 



with the observations. A few studies have considered the effects of turbulent mixing, using an app roach based in mixing-length 
theory, and have suggested that this could also suppress the gas-phase H2O and O2 abundances ( Chieze fc Pineau des Foretsl 
1989; IXie. Allen fc Langeri 1995) However, our current results would appear to rule this out as a solution to the 'water 
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problem'. Another possibility is that there is some as-yet unidentified problem with the reaction rate coefficients used for the 
oxygen gas-phase chemistry, but this is not an issue that our dynamical models can address. 

Since our results for the H2O and O2 abundances are unrealistic compared to those measured in real molecular clouds, 
for one or more of the reasons noted above, will not discuss these molecules any further in this paper, and will focus our 
attention in the following sections on H2 and CO. 



6 DENSITY AND TEMPERATURE PROBABILITY DISTRIBUTION FUNCTIONS 



6.1 Probability density functions (PDFs) at t = tend 

In Figure [7] we plot the mass- weighted and volume- weighted probability density functions of the total gas number density n to t 
at the end of simulations Rl, R2 and R3, i.e. at t cn a = 5.7Myr. We note first that both PDFs have a log-normal shape around 
their peak, although they deviate from this shape in the far wings of the distribution. This form for the density distribution is 
not unexpected. Previous studies o f the density PDF produced by fully-developed supersonic turbulence in isothermal gas find 
that i t has a log- normal form (e.g. Padoan et all 1997: Passot fc Vazquez-Semadeni 1998 : Nordlund fc Padoan | 19991 : Klessen 



2OO0l : lOstriker. Stone fc Gammiell200ll ; IlI et al.ll2004l : iLemaster fe StonelboOsJ : iFederrath. Klessen fc SchmiddboQgh : 



Pa ds 



V2^ 



: exp 



2a] 



ds, 



(26) 



?/2 



where s — ln(p/ (p)), (p) is the mean density of the gas, and where the mean (s) is related to the dispersion a s by (s) = 
due to the constraint of mass conservation. However, the tails of the PDFs may significantly deviate from this log-normal 
dist ribution due to inter mittent fluctuations ( Kritsuk et al. 2007 : Schmidt et al. 2009 : Federrath et al. 20091 ) . 



Padoan et all (1997) argue that the logarithmic density dispersion a 3 is related to the RMS Mach number of the flow, 



M, by 



\n(l + b 2 M 2 ) 



(27) 

where b ~ 0.5. More recently, Federrath. Klessen fc Schmidt ( 20081 ) have shown that the proportionality parameter b depends 
on the relative strength of solenoidal compared to compressive modes in the forcing field used to initialize and drive the 
turbulence. For purely solenoidal forcing, they find that 6 ~ 1/3, while for purely compressive forcing, fo ~ 1. 

We have measured the proportionality constant b using equation (|27[l in the regime of fully developed turbulence 
(3.8 Myr < t < t on( j) and find a mean value of b = 0.32. However, b decreases systematically in time from b = 0.36 at 
t = 3.8 Myr to b = 0.27 at t = t cnd . 

Instead of using the total number density to estimate b, we additionally used the PDF of CO number density to compute 
b. Since the PDF of CO number density significantly departs from a log-normal distribution, we make use of the expression 
for the linear density dispersion I Padoan et al. 1997 : Passot fc Vazquez-Semadeni 19981 : Federrath et al. 20091 ). 



o-p = bM . 



(28) 



This equation for the dispersion does not assume a log-normal distribution ( Federrath. Klessen fc Schmidt 20081 ) . Again, we 
find a systematic decrease of b in time with b = 0.39 at t — 3.8 Myr and b — 0.33 at t — t en d. The mean value in the regime 
of fully developed turbulence is b = 0.35, slightly larger than our estimate using the total number density. This is most likely 
due to the broad plateau of small CO number density seen in its PDF (see Fig. [9}. 

Although the gas in our simulations is non-isothermal, the deviations from isothermality do not ap pear to cause major 
change s in the density PDF compared to the isothermal case. This is consistent with the previous findings of lGlover fc Mac Low 



1 2007bl ) for gas at a slightly lower mean de nsity (no = 100 cm in th e majority of their runs, compared with no = 300 cm 



here). However, simulations such as that of lHennebelle fc Auditl (|2007l ) that consider much lower mean densities and so probe 
the regime in which the gas is thermally unstable produce a broad, bimodal PDF that is not log-normal. 

Regarding the numerical convergence of the density PDF, we note that we find good convergence over much of the 
density range probed by our simulations, but that there is not yet convergence in the wings of the distribution (see also 
Hennebelle fc AuditlE)07l : iFederrath et al.ll2009h . In particular, there appears to be a systematic shift to higher densities with 
increasing numerical resolution that is particularly apparent in the mass- weighted version of the plot, and that causes a slight 
shift in the position of the peak. This behaviour is to be expected, since w e are u nable to fully resolve the shocks in our 
simulations, even at our highest numerical resolution. As lGlover fc Mac Lowl (|2007ah demonstrate, the characteristic cooling 
length of shock-heated gas, L C ooi, at our mean density is of the order of 0.01 pc, and so to resolve these cooling lengths with 
four grid cells in our 20 pc box, we would need to use a numerical resolution of 8000 3 , far larger than is currently possible 
(although [Hennebelle fc Auditll2007l have successfully resolved the post-shock cooling regions in two-dimensional simulations 
of interstellar turbulence). To properly resolve shocks occuring in gas denser than the mean, where the cooling length is 
smaller, we would require an even higher numerical resolution. Since our shocks are under-resolved, the effect of increasing 
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resolution study resolution study 




Figure 7. (a) Mass- weighted PDF of total number density ntot at a time t = i cn( j in runs Rl (green line), R2 (blue line) and R3 (red 
line), (b) As (a), but for the volume- weighted PDF. 
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Figure 8. (a) Mass-weighted PDF of H2 number density tih 2 at a time t = t cn< j in runs Rl (green line), R2 (blue line) and R3 (red 
line), (b) As (a), but for the volume- weighted PDF. 



the numerical resolution is to decrease the size of the cooling regions behind the shocks (which cannot be smaller than the 
grid spacing Ax, even if L coo \ <C Ax), which also allows for greater compression of the cold, post-shock gas. Fortunately, 
the effect of this on the density PDF appears to be relatively small, and hence we should be able to trust the results of 
our simulations, provided that the quantities of interest are not dominated by the behaviour in the wings of the density 
PDF. We also note that quantities that are dependent on the behaviour of the wings of the PDF are in any case difficult to 
characterise based on a single simulation, since the wings a re strongly affected by turbulent intermittency ijKritsuk et al] |2007: 
Federrath. Klessen fc Schmidtj[2009l : iFederrath et al.l [2009 ) , and so vary from realisation to realisation of the same turbulence 



simulation, even if the overall shape of the PDF remains the same when averaged over long enough times. 

In Figure [8] we plot the mass- weighted and volume- weighted PDFs of the H2 number density nn 2 at t = t cn( j. Comparing 
these with the PDFs of the total number density shown in the previous Figure, we see that there is very good agreement. 
This is to be expected: given that only a few percent of the gas remains in atomic form at this point in the simulations, it is 
unsurprising that the PDF of the H2 number density closely follows the PDF of total number density. 

In Figure [9] we plot the mass- weighted and volume- weighted PDFs of the CO number density nco at t = t cn d- Unlike 
the PDFs of total number density and H2 number density, this is not lognormal. The mass-weighted PDF has a clear peak at 
nco ~ 1CF 1 cm -3 , and falls off sharply at higher CO number densities in a fashion similar to a lognormal, but at lower CO 
number densities there is a clear feature in the distribution function, which decreases only slightly from nco ~ cm -3 
down to nco ~ 10 -8 cm -3 . The volume-weighted PDF also shows a peak at nco ~ 10 _1 cm -3 , but is actually bimodal, 
with a second peak at nco ~ 10 -7 cm -3 , although it is also clear that a large number of zones have CO number densities in 
between these two values. 
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Figure 9. (a) Mass-weighted PDF of CO number density nco at a time t = i cnt j in runs Rl (green line), R2 (blue line) and R3 (red 
line), (b) As (a), but for the volume- weighted PDF. 




The curious shape of the CO number density PDFs can be better understood once we realise that we are dealing with 
a PDF made up of two separate contributions, one coming from grid zones with xco — 2;c,tot and a second from grid zones 
with xco £c,tot, and that the dependence of nco on ntot is very different in these two sets of zones. 

If xco — xc,%ot, or in other words if almost all of the carbon in a zone has been converted to CO, then clearly nco 
is directly proportional to ntot- The contribution that these zones make to the PDF therefore simply mirrors the lognormal 
shape of the underlying mass density PDF. A significant fraction of the carbon in our simulation is located in such fully 
molecular zones, and it is the contribution of the gas in these regions that gives us our high density peak in the CO number 
density PDF. 

On the other hand, if xco <C xc,%ot and, crucially, if the CO fraction is not tightly correlated with the total gas density, 
then we would expect to find only a weak correlation between nco and ntot- In other words, if the scatter in the values of 
a^co in gas with a given n to t is large, then this scatter will wash out the effects of any correlation between nco and ntot- As 
we will see later, in Section [Jj a considerable fraction of the gas in our simulations behaves in this fashion, and it is this that 
is responsible for the extended low-density plateau that we see in the CO number density PDF. We will return to this point 
in Section [Jj 

As far as the numerical convergence of the PDFs is concerned, we see very good agreement between the results of our 
three runs over a very wide range of CO number densities. Differences between the three runs are only apparent in the tails 
of the distribution. Increasing the resolution increases the mass fraction in regions with very high CO number density, which 
simply reflects the fact that we better resolve the dense, post-shock gas in the highest resolution simulation, as already noted 
above. We also find fewer regions with very low CO number density in our higher resolution simulations, which again seems 
to be a consequence of the resolution-dependence of the wings of the underlying mass-density PDF. 
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Finally, in Figure [lOl we plot the mass- weighted and volume- weighted PDFs of the gas temperature T at t = t cn d- Several 
features of these plots stand out. First, as the temperature approaches 10 K, the PDF falls steeply, and there is almost no gas 
in the simulation with T < 10 K. This feature of the plot is artificial, and is a consequence of our adoption of a temperature 
floor at T = Td us t = 10 K in our treatment of the radiative cooling (see £12.3.21 above) . Thus, the only fluid elements that 
have T < 10 K are those that had temperatures close to 10 K and then underwent a strong rarefaction, leading to significant 
adiabatic cooling. Moreover, this must have happened within the last 0.1 Myr, or else cosmic ray heating would have warmed 
the gas up above 10 K again. As is apparent from the PDFs, very few of the fluid elements in our simulations find themselves 
in this situation at any given time. 

The second obvious feature of the temperature PDFs is the clear power-law tail between 30 K and 120 K in both the 
mass and the volume-weighted PDFs. This tail is composed of gas with a low Ay >e g that is heated primarily by photoelectric 
emission from dust grains and is cooled primarily by C + fine stru cture emission. The equilibrium tempera ture of this gas 
varies approximately as a power-law function of density, T cq oc n ' 7 ( Larsonll2005l ; iGlover fc Mac Lowll2007bl ). 

Regarding numerical convergence, we again find good convergence for the majority of the PDF, with significant differences 
visible only around the low temperature peak in the distribution. Increasing the numerical resolution shifts the peak to slightly 
lower gas temperatures, reflecting the fact that the coldest gas is also, typically, the densest, and that this dense material is 
better resolved in the higher resolution simulations. 



6.2 Time evolution of the PDFs 

In Figure [TTk . we show how the mass- weighted PDF of the total number density evolves with time in the 256 run R3. The 
first output time for which data is plotted, t — 0.6 Myr, corresponds to less than half of a turbulent crossing time, and at this 
point in the simulation, the imprint of the initial conditions is still quite apparent. The density PDF at this time is clearly 
not lognormal. Instead, it is bimodal, with one peak close to the starting density no = 300 cm -3 , corresponding to gas which 
has not yet been significantly affected by the turbulence, and a second peak at n ~ 3000 cm -3 , corresponding to gas that has 
already been compressed by the strong, large-scale shocks present in the initial turbulent velocity field. 

By the time of the second output dump at t — 1.9 Myr, corresponding to roughly one turbulent crossing time, the picture 
is quite different. The characteristic lognormal density PDF has now been established, although a few fluctuations in the low 
density tail of the PDF are still apparent. These have vanished by the time of our third output dump, at t = 3.2 Myr, and 
from this point on we see very little evidence for any change in the PDF, suggesting that the density distribution has reached 
a statistical steady state. 

In Figure lllb . we plot the mass-weighted PDF of the H2 number density in run R3 at various output times. As in 
Figure [TTk . the PDF at the earliest output time is clearly not lognormal, but after one turbulent crossing time, it has become 
lognormal over much of the range of densities plotted. This is easily understood when one considers that by t — 1.9Myr, roughly 
80% of the hydrogen in the simulation has already become molecular (see Figu re [4j and that we expect there to be a clear 



correlation between gas density and molecular fraction (see Glover fc Mac Low 2007bl . or jJ7] below). At gas densities where 



the hydrogen is almost fully molecular, we expect that the H2 number density PDF will simply track the underlying density 
PDF. Comparison of Figures [TTk and lllb shows that this is the case for H2 number densities greater than riH 2 ~ 10 cm -3 
at t = 1.9 Myr. Lower H2 number densities correspond to regions of low density gas where a significant fraction of atomic 
hydrogen remains, and in these regions there is no simple mapping between gas density and H2 number density, since here 
the H2 fraction is sensitive not only to the current density and th e degree of shielding, but also to the previous dynamical 
history of the gas ( Glover fc Mac Low 2007b ; Federrath et al. 2008]). It is therefore not surprising that at these low densities 



the PDF is not lognormal. 

At later output times, we see little change in the PDF at densities nu 2 > 10 cm -3 , which simply reflects the fact that 
there is very little change in the underlying density PDF. At lower densities, the PDF continues to evolve with time. The 
mass fraction of gas with extremely low H2 number densities, tih 2 < 0.3 cm -3 , significantly decreases. This gas physically 
corresponds to gas with a low total number density and low H2 fraction. As time passes, the H2 fraction in this gas increases, 
and so the H2 number densities increase, even though the statistical distribution of the total number densities is unchanged. 
At the latest output time, t = 5.7 Myr, the shape of the full H2 number density PDF is almost identical to that of the full 
density PDF, since by this point, almost all of the hydrogen has become molecular. 

In Figure [TTb . we show how the mass- weighted PDF of the CO number density evolves with time in run R3. Just as with 
n to t and iih 2 , the shape of the PDF at our earliest output time bears little resemblance to its later shape, while by the time 
of the second output dump, the shape of the final PDF has become far better established. Unlike ntot and nn 2 , however, we 
continue to see evolution in both the low and the high density portions of the PDF, and it is only after 3.2 Myr, corresponding 
to 1.6 turbulent crossing times, that the CO number density PDF reaches a statistical steady state. This longer evolutionary 
timescale is a consequence of the longer time required to form CO, compared to H2. Between t = 1.9 Myr and t = 3.2 Myr, 
the total mass of H2 in the simulation increases by about 15%, while the total mass of CO increases by about 50%. 

Finally, in Figure [TTH . we show how the mass- weighted temperature PDF varies with time. At the earliest output time, 
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Figure 11. (a) Evolution with time of the mass- weighted PDF of the total number density ntot in our 256 3 zone run, R3. Values are 
plotted for t = 0.6Myr (green long-dashed line), t = 1.9Myr (black dotted line), t = 3.2Myr (blue short-dashed line), t = 4.4Myr (mauve 
dot-dashed line) and t = 5.7 Myr (red solid line), (b) As (a), but for the H2 number density, (c) As (a), but for the CO number density, 
(d) As (a), but for the gas temperature. 



this has a two-peaked structure: a high temperature peak, centered on T ~ 65 K, corresponding to unshocked gas that still 
has a density and temperature close to its initial value, and a low temperature peak, centered on T ~ 30 K, corresponding to 
shocked, higher-density gas (compare with the density PDF at t = 0.6Myr in Figure [TTV ). By t = 1.9 Myr, this double-peaked 
structure has disappeared, but as with the CO number density PDF, it is not until t = 3.2 Myr that large portions of the PDF 
settle into their final form. By this point, the power-law high temperature tail has become fully established, and the PDF at 
temperatures T > 30 K shows very little evolution at later times. On the other hand, the portion of the PDF around the low 
temperature peak continues to evolve. The PDF peaks at T ~ 16 K at t = 3.2 Myr, but by t — 4.4 Myr the peak has shifted 
to T ~ 13 K, while by t — 5.7 Myr, it has shifted further, to T ~ 12 K. We therefore see a systematic increase with time of the 
amount of cold gas present in the simulation. This increase has not yet come to an end by the end of our simulation, although 
it has significantly slowed: there is far more evolution in the temperature PDF between t = 3.2 Myr and t = 4.4 Myr than 
there is between t = 4.4 Myr and t = 5.7 Myr. This increase in the mass of cold gas is simply driven by radiative cooling of 
dense, fully-molecular gas. At temperatures T> 10 K, the cooling time of the gas is very short, but for temperatures close to 



our t emperature floor of 10 K, it becomes comparable to the dynamical time (see e.g. Figure 2a in lNeufeld. Lepp fc Melnick 
1995). 



7 SPATIAL DISTRIBUTION OF H 2 AND CO 

In Figure [12k . we plot the column density of hydrogen nuclei, A^H.tot, projected along the z-axis of the simulation (i.e. along 
the axis parallel to the direction of the initial magnetic field) in run R3 and at time t = t cn d- In Figures [12b and 112b . we 
show similar plots of the H2 and CO column densities. The basic morphology of the gas is the same in all three plots. The 
gas has a filamentary distribution, and large spatial variations in the column densities are apparent, including, coincidentally, 
a rather prominent under-density visible toward the bottom-left of the figures. The plots of A^H.tot and Nh 2 are very difficult 
to distinguish, which is unsurprising since most of the hydrogen gas in the simulation is molecular by this point. On the other 
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hand, the plot of CO column density is clearly different from the other two plots: the underdense regions are larger, and also 
more numerous, particularly towards the edges of the box. Similar results are found if we consider projections along lines of 
sight perpendicular to the initial magnetic field. 

The difference between the spatial distributions of Nco and iVn 2 can be more clearly highlighted by plotting the ratio 
Nh 2 /Nco , as we do in Figure 1121 1. Along lines of sight corresponding to the highest column densities, this ratio is around 
10 3 ' 5 , as we would expect for gas in which all of the hydrogen is in the form of H2 and all of the carbon is in the form of 
CO. However, if we look along lines of sight that pass through regions of lower total column density, then we find values for 
this ratio that are up to a factor of thirty larger. Along these lines of sight, much of the carbon remains in the form of C 
or C + . It is clear from Figure [I2h that the regions with low CO column densities are found preferentially towards the edges 
of the simulation volume. This is only to be expected, given our treatment of the external ultraviolet radiation. Gas close 
to an edge of the simulation volume will tend to have a low column density of dust between itself and the edge, and so will 
be more readily affected by UV photons propagating inwards in that direction, even if in other directions, that same gas 
is well-shielded. However, it should be stressed that this is not simply an edge effect: when the column density is low, UV 
photons can penetrate in to the volume to considerable depths, and are not confined to a narrow region at the surface of the 
box. 

Similarly, the two large regions with high H2:CO column density ratios that are close to the center of the projections are 
well-shielded by dust along the x and y axes of the simulation, but have a low column density of material shielding them in 
the z direction, and so are strongly affected by UV photons propagating inwards in that direction. 

Clearly, if we were to solve for the UV radiation field using much higher angular resolution, we would expect to obtain 
slightly different results for the spatial distributions of the CO column density and the H2:CO column density ratio. Never- 
theless, the picture we obtain here is qualitatively correct. The dumpiness created by the turbulence opens up channels in 
the gas distribution, allowing UV photons to propagate far deeper into the cloud than would be possible if it were a single 
homogeneous mass of gas (see also Boisse 199pl : Padoan et al. 20041 : Bethell. Zweibel fc Li 2007 '). 



7.1 CO to H 2 ratio 

The results shown in Figure [12] indicate that the amount of CO in the gas depends in part upon the column density of the 
material shielding the gas, but being two-dimensional projections, they do not easily allow us to quantify this dependence, 
or to separate the effects of increased shielding from the effects of increased gas number density. To do this, we need to make 
use of the full three-dimensional spatial information contained in the datacube. 

In Figure [13b . we show how the fractional abundance of CO, xco, depends on the number density n to t by plotting the 
mass- weighted two-dimensional PDF of these two quantities. In Figure [13b . we show a similar plot of xco versus the effective 
visual extinction Av,oft . For a grid cell at position x, we define the effective visual extinction Av.eff as 



-2.5A v (x p+ ) , -2.5A v (i p _) 



: e ~ vi.-.p-w + e 



(29) 



where Ay(x p +) is the visual extinction of material between the cell and the edge of the volume in the positive direction 
along the x p axis, and Ay{x p -) is the same in the negative direction. The choice of the factor of 2.5 occurs because the CO 
photodissociation rate scales with the visual extinction Ay as exp(— 2.5Ay). The value of Ay, off corresponds to the visual 
extinction used in our code, in the context of our six-ray approximation, for computing the CO photodissociation rate (see 
Section |2~2[) . 

Figure [13k demonstrates that there is only a weak correlation between the CO abundance and the number density. At 
densities ?i t ot < 10 cm" 3 , most of the gas has a very low CO abundance, xco < 10~ 8 . On the other hand, at densities 
ntot > 10 4 cm~ 3 , almost all of the carbon in the gas is found in the form of CO. At densities between these two extremes, 
however, there is a very wide spread in :rco at any given n to t. For instance, at n to t = 100 cm -3 , values of xco between 10 -9 
and 1.4 x 10~ 4 are almost equally probable. This large scatter in xco is unexpected and demands an explanation. 

A large part of the puzzle is explained by Figure [13b . which shows that there is a correlation between xco and the effective 
visual extinction, Ay >c g, that is stronger than the correlation between CO fraction and gas number density. This is in line with 
our expectations based on our discussion of the H2:CO column density ratio above, and is also readily explained by considering 
the microphysics of CO formation and destruction. Consider the following simplified model for the CO abundance. CO is 
formed from gas phase C or C + by a variety of reactions, but the most important ones involve either a hydrocarbon radical 
(e.g. CH) or its ion (e.g. CH + ), or the OH radical or its ion, OH + . In cold gas, the formation of one of these intermediate 
species is the rate-limiting step in forming CO, as the various routes by which these species can be formed typically involve 
either a radiative association reaction with H2, or the cosmic ray ionization of H2, both of which are slow processes. Let us 
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Figure 12. (a) Column density of hydrogen nuclei, A'jj tot , in run R3 at time t = t cn< i, viewed along a line of sight parallel to the z-axis 
of the simulation volume. This direction is also parallel to the initial orientation of the magnetic field, (b) As (a), but for the H2 column 
density, (c) As (a), but for the CO column density, (d) Ratio of H2 column density to CO column density along the same line of sight 
through run R3 at time t = t en( j . 



suppose, for simplicity, that reactions involving hydrocarbon radicals and ions dominated In that case, we can write the CO 
formation rate as -Rform"-c.totn.H 2 , where nc.tot = nc + n c + , and where R{ olnl is the formation rate of our intermediate species, 
multiplied by a factor that accounts for the fact that some of the intermediate radicals and ions will be photodissociated, 
rather than reacting to form CO (or a further intermediate, such as CO + , that reacts rapidly to form CO). If CO is primarily 
destroyed by photodissociation, at a rate RpdTico, then in chemical equilibrium, the CO fractional abundance is given by 



-Rforr 

~R~ 



P d 



(30) 



The photodissociation rate R p d can be written in terms of Ay.eff as R p d = 2 x 10 _10 / 8 h exp(— 2.5Av, c ff ), where / s h = /co/h 2 
is the product of the shielding factors due to CO self-shielding (fco) an d due to the shielding of CO by H2 (/h 2 ) that we 
introduced in H2.2I We can therefore rewrite Equation 1301 as 

Rio 



xco = 



( Rtor m \ r - 

\2x 10- 10 ) /sl 



£C,totnH 2 - 



,2 x 10- l! 

Consideration of the different processes involving C or C + that lead to the formation of CH or CH + suggests that R{ 



(31) 



We could construct a very similar model in the case that reactions with OH and OH + dominate, only with the number density of 
atomic oxygen, no, playing the role of np tot above. 



20 Glover et al. 



t = 5.7 Myr, mass-weighted PDF log 10 PDF t = 5.7 Myr, mass-weighted PDF log 10 PDF 




O 
o 



o) -9 



-10 



t = 5.7 Myr, mass-weighted PDF 



log 10 PDF 

n 

-1 

-2 
-3 
-4 
-5 
-6 



t = 5.7 Myr, mass-weighted PDF 



O 
o 



o) -9 



-10 



12 3 4 
lo 9io n tot exp[+2.5 A V etf ] 



2 3 4 5 

'ogirjntotW 1 exp[+2.5 A Vieff ] 



log 10 PDF 

rn 

-1 
-2 
-3 
-4 
-5 
-6 



Figure 13. (a) Mass- weighted two-dimensional PDF of CO fractional abundance ico versus gas number density ntot- (b) Mass- 
weighted two-dimensional PDF of CO fractional abundance ico versus effective visual extinction A\r c ft (defined in Equation I29|l . (c) 
Mass-weighted two-dimensional PDF of xco versus exp(+2.5j4v O ff)ntot- (d) Mass-weighted two-dimensional PDF of ico versus §, 
defined in Equation 1341 



should have a value of roughly 10 17 cm 3 s 1 , give or take an order of magnitude. If we assume, again for simplicity, that in 
fact iiform = 2 x 1CF 17 cm 3 s -1 , then Equation 1311 becomes 

xco = 10-V- 5Av .° ff / s 7>c,totn H2 . (32) 

If we further simplify matters by assuming that the H2 fraction is of order one, and that most of the carbon in the gas is still 
in the form of C or C + , so that a;c,tot ~ 10 -4 , then we obtain the following expression for xco 

xco ~ lO^V-^-^/^n. (33) 

Given the large number of assumptions and simplifications that we have made above, this expression should clearly be 
treated only as a rough order-of-magnitude estimate of the CO abundance. Nevertheless, this simplified model does highlight 
some of the behaviour that we see for the actual CO abundance in our simulations. Our simplified model predicts that the 
abundance of CO should vary only linearly with the gas number density, and with /7 (which itself is a complicated function 
of Nn 2 and Nco), but should vary exponentially with ^v.cff • Thus, small changes in Ay, c ff will produce a much larger change 
in xco than even relatively large changes in n, and so we would expect to find a much stronger correlation between xco and 
^4v,eff than between xco and n to t, as indeed is the case in our simulations (see Figure [13]). 

If the behaviour of xco is determined primarily by the amount of shielding by dust and molecules, rather than by the 
density, then we would expect to obtain a much tighter correlation if, instead of plotting ico against n to t , we plot xco against 
ntot exp(2.5v4v,eff ) or, even better, against 

£ = n tot exp(2.5^4v, e ff)/~ h 1 . (34) 

This we have done in Figures 113b and I13H . These Figures show that when the CO fraction is small, then we obtain a much 
tighter correlation between xco and £ than between xco and n to t, while accounting for the effects of shielding from CO and 
H2 improves the correlation even further. This confirms that most of the scatter we see in the a;co _ ntot plot at low xco is 
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Figure 14. Mass- weighted two-dimensional PDF of iy.cff versus ntot. 



due to the fact that the CO fraction is far more sensitive to changes in the amount of shielding than to changes in the density, 
and that the amount of shielding is not well correlated with the gas density. This point is further emphasised in Figure 1141 
which shows that there is indeed a large scatter in the effective visual extinction at most densities. 

If the CO abundance is large, then the relationship between xco and the amount of shielding becomes quite different, and 
the clear correlation found at low xco rapidly vanishes. Our simplified model for xco again can give us insight into why this 
happens. Consider that in order to produce a CO abundance xco ~ 1CP 5 , our model predicts that we must have £ ~ 10 6 . For 
a gas density n ~ 300 cm -3 , this means that the CO photodissociation rate is reduced by a factor of roughly 3000 compared 
to the optically thin value, i.e. 7? p d ~ 10~ 13 s -1 . The corresponding photodissociation timescale is then 0.3 Myr, which is not 
negligible in comparison to the dynamical timescale of the gas, and so it is probably no longer safe to assume either that the 
gas is in chemical equilibrium or that photodissociation dominates the destruction of CO. As our model was based on both 
of these assumptions, it is unsurprising that it breaks down at high xco- 

To summarise, what these results are showing us is that in low extinction gas, the CO abundance is determined by a 
combination of the extinction and the density, with the extinction playing the primary role, while at high extinction, other 
physics, such as the dynamical history of the gas, or the impact of other destruction mechanisms, such as charge transfer with 
He + , becomes far more important. 

The relationship between Av.eff and ntot, shown in Figure [141 is also worth discussing. As noted above, there is no strong 
correlation between these two quantities. However, they are not completely uncorrelated either. There is a clear deficit of 
zones with high -Av.eff at low densities, and there is also a deficit of zones with low Ay. c s at high densities. The latter is a 
consequence of the finite resolution of our simulation, since the value of Av.efi in a given zone has a lower bound 

0.5nAa; 

Av - cS - 1.87 x lO^i cm-2 ' (35) 
where Ax is the size of a grid zone, corresponding to absorption within the zone itself. However, this does not explain the 
deficit of points with high extinction and low density. Instead, this is due to the fact that there are only very few voids in the 
density distribution that are entirely surrounded by high extinction material. In addition, little of the mass in the simulation 
volume is to be found at these very low densities. 

Also of interest is the fact that the smallest values of Ay, e ff found in our simulation are roughly 0.1-0.3, rather than zero. 
The reason for this is quite simple. The regions in the simulation with the lowest values of Av, c s are those at the edges and 
corners. These regions are highly exposed to radiation propagating inwards from the nearest edges of the box. However, at 
the same time, they are very well shielded from radiation propagating inwards from the opposite edges of the box, since they 
are protected by the full width of gas in the simulation. Therefore, these zones see a mean intensity that is smaller than the 
mean intensity that they would see if the opacity in every direction were zero, and hence they have a non-zero Av, e ff, even 
though the visual extinction between them and the closest edge of the box may be very close to zero. It is also clear that 
our "six-ray" approximation tends to overestimate the degree to which the zones at the edges of the box are shielded. For 
instance, consider a zone located at the edge of the box, in the center of one of the six faces. In one direction it has a very low 
Ay, however, in the other five directions, it may have a very high Ay. Therefore, using our approximation we would predict 
that it would see a mean intensity of radiation that was roughly one-sixth of the value in the absence of shielding, while in 
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Figure 15. Mass-weighted two-dimensional PDF of the fractional abundance of atomic carbon, xq versus effective visual extinction 

Av.eff ■ 



reality it should see a mean intensity that is only one-half of the optically thin value. Fortunately, this is only really a problem 
for zones right at the edge of our simulation, and most of the gas should not be significantly affected. To test this, we verified 
that the removal of 32 zones of material from each of the edges of our 256 3 simulation did not significantly affect the PDFs 
reported in this paper. 



7.2 Comparison to photodissociation region models 



To conclude this section, it is interesting to compare our results for the CO distribution with what we would expect based 
on classical models of uniform density photodissociation regions (see e.g. iHollenbach fc Tielensiri999l . and references therein). 
These predict that we should find an outer shell of C + surrounding a shell of neutral carbon, that in turn surrounds a central 
region dominated by CO. If we examine how the C and CO abundances vary with visual extinction (Figures ll3b and !15|l . then 
they are to some degree consistent with this picture. We see a peak in the typical atomic carbon abundance at Av, e fl ~ 1, 
followed by a decline at higher .Av.cff ■, as CO begins to dominate. However, what is apparent in both figures is the large degree 
of scatter around this underlying trend. There are many regions of high extinction in which the C abundance is very small, as 
the PDR models would predict, but at the same time there are also regions in which xc remains large at high Ay. Similarly, 
although most regions with Ay > 2 are CO dominated, there is also a significant amount of mass with high extinction that 
still has a very low CO abundance. 

Some of this scatter is surely due to the highly inhomogeneous nature of the density field. As we have already seen, the 
density of gas with a given ^4y,off can vary over several orders of magnitude. In addition, it is likely that turbulent mixing 
of material between high extinction and low extinction sightlines also plays an impo rtant role in determining the spatial 
distribution of C + , C and CO (see also Glover fc Mac Low 2007bl : Federrath et al. 200ct ). To properly disentangle these effects 
lies beyond the scope of this introductory paper, but is something we plan to revisit in future work. 



8 SUMMARY 

In this paper, we have outlined the methods that we have developed to model the coupled thermal, chemical and dynamical 
evolution of the turbulent gas making up giant molecular clouds. We have shown that it is now possible to perform high- 
resolution simulations in three dimensions that track the thermal and chemical evolution of the gas in a realistic fashion. The 
chemical network that we use in our simulations is significantly simplified compared to the most detailed models available, 
but nevertheless performs with acceptable accuracy for our purposes, and the largest errors in our models come not from our 
simplified chemistry, but from our approximate treatment of the effects of the external radiation field. 

We have performed simulations with numerical resolutions of 64 3 , 128 3 and 256 3 , and have demonstrated that most of 
our results are well-converged in our 256 3 run. The main exception is the high density tail of the density PDF, which is not 
fully converged at a resolution of 256 3 zones (although the peak of the density PDF is well converged). This results from 
the improvement of our ability to resolve dense post-shock gas with increasing numerical resolution, and is also responsible 
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for the small dependence on resolution that remains in our determination of the O2 and H2O abundances, and in the gas 
temperature distribution at the lowest temperatures. 

We have also quantified the timescales on which various quantities of interest reach a statistical steady state. The density 
PDF, which has the same log-normal form as found in many previous simulations of interstellar turbulence, reaches a steady 
state after roughly one and a half turbulent crossing times, as does the H2 number density distribution in all but the lowest 
density gas. On the other hand, the CO number density distribution takes closer to two crossing times to settle into a steady 
state, and the temperature distribution still has not reached a final steady state after three crossing times. 

We find that CO formation occurs rapidly in the dense, turbulent gas studied here. Most of the CO is produced within 
the first 2-3 Myr of the simulation, i.e. within 1-2 turbulent crossing times. For comparison, most of the H2 produced in 
the simulation forms within a single crossing time. These short chemical timescales suggest that the limiting timescale in the 
formation of a molecular cloud is not the time required to convert the hydrogen to H2 and the carbon to CO, but rather the 
time required to assemble the cloud material from the low density ISM. Once large enough spatial and column densities are 
reached, conversion of the gas to molecular form should follow very quickly. 

Our simulations also demonstrate that the gas in molecular clouds is not isothermal. Most of the molecular gas has a low 
temperature, in the range of 10-20 K, but there is also a clear power-law tail in the temperature PDF at higher temperatures, 
made up of low extinction gas heated by photoelectric emission from dust grains. 

We have found that the CO abundances produced in our simulations are not well correlated with the gas density. This 
means that the use of a density cut to identify molecular regions in turbulence simulations that do not self-consistently follow 
the chemical evolution of the gas is dangerous, and may give misleading results. Similarly, the assumption that all observed 
CO is at roughly the same gas density may be similarly misleading, as is the assumption that all of the CO in a molecular 
cloud is located at densities that are high enough for its rotational lines to be thermalized. 

The poor correlation that we find between CO abundance and gas density is explained by the sensitivity of the CO to 
the amount of UV shielding provided by the dust, H2 and CO. Variations in the amount of shielding produce a far greater 
variation in the CO abundance than do variations in the gas density. In addition, the amount of shielding provided by the 
gas is not well correlated with the density, owing to the highly inhomogeneous nature of the density field. The important role 
tha t density inhomogeneities play in determining the chemical structure of molecular clouds ha s long been recognized (see 
e.g. Istutzki fc GuesteJl990l : Icierens. Stutzki fc Winnewisser|[l992l : I Williams. Blitz fc Starkll995h , but the most sophisticated 
current PDR models continue to rely on a somewhat artificial picture of their structure, visu alising them as an ensemble of 
clum ps with well-defined densities surrounded by a much lower density interclump medium (e.g. lSun et al I I2OO8I ;lKramer et al 



2008). The picture suggested by our simulations is rather different. We find a continuous density distribution, with no clear 



separation between 'clump' and 'interclump' material, and with continual turbulent mixing of gas between low-density and 
high-density regions. We look forward to exploring the observational consequences of this picture in greater detail in future 
work. 
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APPENDIX A: ENSURING CONSISTENT ADVECTION OF CHEMICAL SPECIES 
Al Consistent Multi-fluid Advection (CMA) 

Let Xm(i,j, k) be the fractional abundance relative to hydrogen of chemical species m in grid zone k), let X n (i,j, k) be 
the elemental abundance, relative to hydrogen, of element n = H, He, C or O in grid zone (i, j, k), and let N n ,m be the number 
of nuclei of element n in species m. Then for each element n, we can construct a constraint equation of the form 

^ N ntm x m {i, j, k) = X n {i,j,k). (Al) 

m 

Now, if we assume that the elemental composition of the gas remains the same throughout our simulation volume, i.e. that 
the elemental abundances X„ are independent of positionp then these constraint equations become 

^ N n , m x m (i, j, k) = X„. (A2) 

m 

Although it is trivial to ensure that these constraints are satisfied at the start of our simulations, ensuring that they remain 
satisfied as the gas is evolved forward in time is not so simple. 

The root cause of the problem is to be found in the advection scheme. Most Eulerian treatments of multi-species flows 
model the species abundances x m as passive scalars that are advected with the flow. If we consider a one- dimensional scheme, 
for simplicity, then we can write the advection equation for a given species m in grid zone i in finite volume form as 

piXm,i(t + Ai)AVi = pix m ,i(t)AVi - At [A i+1/2 F i+1/2 - A i ^ 1/2 F i ^ 1/2 ] , (A3) 

where F i+1/2 = Pi+i/ 2 Vi+i/ 2 x m , i+1 / 2 is the flux of x m from zone i into zone i + 1, i ? i _i /2 = Pi_i/2«i-i/2Zm,i-i/2 is the flux 
of Xm from zone i — 1 into zone i, A i _ 1 / 2 and A i+1 / 2 are the areas of the interfaces between zones i — 1 and i, and between 
zones i and i + 1, respectively, and AVi is the zone volume. 

If the elemental abundances X n are independent of position, then at both interfaces of zone i, the species fluxes should 
satisfy the constraint equations 



Pi±l/2«i±l/2 = XnPi±l/2Vi±l/2- (A4) 



2 This is a reasonable assum ption in small-scale simulation s of the local interstellar medium, where we can treat all of the gas as having 
the same metallicity (see e.g. iMever, Jura, fc C ardclli 199g), but will not hold on scales large enough that Galactic metallicity gradients 
become important. Note also that we are only assuming a homogeneous metallicity, not homogeneous chemical abundances. 
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In order to ensure stability, fluxes are often computed from the distribution of the flow variables upstre am of the interface 
( "upwinding" ) using some interpolation scheme. For instance, in ZEUS-MP, which uses Ivan Leer! (|1977l ) advection, the up- 
winde d interpolated value q* of a zone-centered scalar q at the negative interface of zone i is given by (jStone fc Norman 
1992a) 



<7i = 



qi-i + (Axi-i — v i _ 1 / 2 At)(dq i -i/2) 
qi + (An + v t _ 1/2 At)(dqi/2) 



Vi-l/2 > 
u<_X/2 < 



(A5) 



where g; is the zone-centered value of q in zone i, Axi is the size of zone i, and dqi is the monotonized van Leer slope, given 
by 



f 2(A gi _ 1/2 A 9i+1/2 ) 

dqi = < A 9i-i/2+A 9l+ i /2 




Aq i+1 /2Aq i _ 1/2 > 
Aq i+1 / 2 Aq i _ 1/ 2 ^ 



(A6) 



where Aq i+1 / 2 = (g«+i — qi)/Axi. (Note that in Equation IA5I we have neglected any motion of the underlying coordinate 
grid, for simplicit y). Given q* , the flux then fo llows from -Fj_i/2 = Q_* v i-i/2- Also in common usage is the Piecewise Parabolic 
Method (PPM) of lColella fc Woodward! |l98J), which uses higher-order interpolation, and so produces a better quality solution 
(at the cost of increased computational effort). 

Unfortunately, as several authors have noted ( Frvxell. Miiller fc Arnett 19891 : Larrouturoul 1991 ; Plewa fc Miiller 19991 ) , 
the fact that the interpolation profiles in these schemes are constructed independently for each chemical species means that 
there is no guarantee that the resulting fluxes actually satisfy Eq. IA4I In general, they do not do so. Moreover, this violation 
occurs even when the underlying advection scheme is conservative. Thus, even if our initial chemical abundances satisfy the 
constraints represented by Equation IA2I as soon as we begin advecting them with a standard higher-order scheme, these 
constraints will be violated. 

This problem is often deal t with by a renormaliz ation of species abundances following the advection step to ensure that 
Eq, IA2l is satisfied. However, as IPlewa fc Muflerl (|l999h note, this procedure lacks any formal justification and can lead to large 
systematic error s in th e abundances of the least abundant species. It can also destroy the conservative nature of the scheme. 
Plewa fc MMerl dl999l ) suggest an improved way of dealing with this problem, which they term the Consistent Multi-fluid 
Advection (CMA) method. They consider the case in which one has only a single constraint equation 



/]Nl, m Xm(i,j, k ) = 1) 



(A7) 



where we have chosen to set X n — 1 for simplicity. The example they consider is a fluid consisting of many different elements 
(but no composite molecules), whose mass fractions must sum to unity. An alternative example, more in keeping with our 
discussion here, is a gas consisting of only one element (e.g. hydrogen) that can form several different stable chemical species 
(e.g. H + , H~, H, H2). They show that in this case, one can satisfy Equation IA7I while preserving the conservative nature of 
the advection scheme by replacing the interpolated abundances x i ±±/2.m used to construct the fluxes in Equation IA3I with 
the modified abundances 

x i±l/2,m 



Xi±l/2,m — 



^.#£±1/2, m 



(A8) 



In other words, instead of ensuring that the constraint equation is satisfied by normalizing the abundances after the advection 
step, we ensure that it is satisfied by normalizing the fluxes during the advection step. 

When dealing with the PPM code, additional modificatio ns are necessary to av oid problems related to PPM's contact 
discontinuity detection algorithm, as discussed in Section 2.3 of lPlewa fc Miiller (1999). However, these problems do not occur 
with the simpler van Leer advection scheme used in ZEUS-MP, and so we do not discuss them further here. 



A2 Extending the CMA algorithm 



Although IPlewa fc Mullen 's CMA algorithm is simple to implement and is highly effective in practice, its applicability to 
chemically reacting flow is limited. The problem comes from the assumption that the chemical abundances must satisfy only 
a single constraint equation. In a flow containing compounds of multiple elements, such as CO, this is not the case. In such a 
flow, our modified interface abundances for carbon and oxygen must satisfy both 

^N c , m x m = X c (A9) 

m 

and 

'Y J No, m x m =Xo (A10) 
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If we rescale the abundances of those species containing carbon according to 
X c 

Xi±l/2.m = 7T £j±l/2.m (All) 

and those species containing oxygen according to 

Xi±l/2.m = ICr 35»±l/2,m (A-12) 

then we will find, in general, that for species containing both carbon and oxygen, x'i±i/2 m 7^ Xi±i/2,m- We therefore cannot 
use the simple CMA prescription to rescale the interface abundances (and the fluxes) in this kind of flow. 

To avoid this problem, we propose an extension of the CMA algorithm, which we term modified CMA (or MCMA, for 
short). The basic idea behind this algorithm is the same as that motivating the CMA algorithm: we aim to modify the chemical 
abundances used to construct the species fluxes, such that the fluxes satisfy the appropriate chemical abundance constraint 
equations, while the advection scheme itself remains conservative. We accomplish this by writing the rescaled abundances as 

Xi±l/2,m = ( ^ Vk I X i± y 2>m (A13) 

V /c " tot . m / 

where Nk, m is the number of nuclei of element k in species m, N to t,m is the total number of nuclei in species m, and the rjk 
are correction factors chosen so that the set of constraint equations 

^2 N k,mXi±l/2,m = X h (A14) 
m 

are simultaneously satisfied. The required correction factors rjt can be found by solving the matrix equation 

Mr; = X (A15) 

where T] = (771 , T]k) , X = (Xi , X^ ) , and where the elements of the matrix M are given by 

M u = V ^h^ NhmXt±1/2tmj (A16) 

* ' ivtot.m 

m 

where we sum over all species m. In a flow with a mix of species containing only a single element, iVi, m /iVtot,m = 1 for all m, 



and this prescription reduces to lPlewa fc Mullen 's CMA scheme 



The main drawback of this scheme compared to the original iPlewa fc Miillerl (1999) scheme is that we must now perform 
a matrix inversion for every flux on every advection step. However, the rank of this matrix scales only as the number of 
elements, not the number of species, and so this additional cost will generally be dwarfed by the cost of solving the chemical 
rate equations themselves. 
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Table Bl. List of collisional gas-phase reactions included in our chemical model 



No. Reaction 



Rate coefficient (cm 3 s 1 ) 



Notes 



Ref. 



1 H + c- -> H" +7 

2 H-+H^H 2 +c- 

3 H + H+ -> H+ + 7 

4 H + H+ -» H 2 + H+ 

5 H-+H+^H + H 

6 H++c-^H + H 

7 H 2 + H+ ->• H+ + H 



8 H 2 + e- ^H + H + c- 

9 H 2 +H^H + H + H 



10 H 2 + H 2 ->■ H 2 + H + H 



11 H + c~ -> H+ +< 



12 H+ + e~ ^H + 7 



13 



H + c" 



fcl 

fe 2 

k:i 

k4 
^5 

fc 6 

fc 7 



: dex[-17.845 + 0.762 logT + 0.1523(logT) 2 

-0.03274(log T) :i ] 
: dcx[-16.420 + 0.1998(log T) 2 

- 5.447 x 10- 3 (logT) 4 

+ 4.0415 x 10- 5 (logT) 6 ] 

: 1.5 X 10~ 9 

: 4.0 x io- 9 T-° 17 
: dcx[-19.38 - 1.523 logT 
1.118(logT) 2 - 0.1269(logT) 3 ] 
: 6.4 x 10" 10 

: 2.4 x 10" 6 T- 1 /2(i.o + T/20000) 
: 1.0 x 10~ 8 
: 1.32 x io- 6 T-°- 76 
: [-3.3232183 x 10~ 7 
+ 3.3735382 x 10~ 7 lnT 

- 1.4491368 x 10- 7 (lnT) 2 
+ 3.4172805 x 10- 8 (lnT) 3 
-4.7813720 x 10~ 9 (lnT) 4 
+ 3.9731542 x 10" 10 (lnT) 5 

- 1.8171411 x 10 _11 (lnT) 6 
+ 3.5311932 x 10- 13 (lnT) 7 ] 
x cxp ( - 212 T 3715 ) 

: 3.73 x 1Q-9T" 1121 cxp (- 



T < 6000 K 



T > 6000 K 
T ^ 300 K 
T > 300 K 



T ^ 617 K 
T > 617 K 



fc 8 

fc 9 ,i = 6.67 x 10" 12 T 1/2 
k a h = 3.52 x 10~ 9 



*9,h 
«cr,H = dex 



cxp | 



cxp 
43900 



[_ ( l + SSS8fi)] 



3.0 - 0.416 log 

996xlO~' ,u T 



.10000/ 
30^4.1881 / 

rasBrexp^- 



0.327 {log (j^)}' 



(1.0+6.761XlO- b T) 

fcio, h = 1.3xl0- 9 cxp (-SMpo) 

™cr,H 2 = dex 
fell 



54657.4 > 
T , 



4.845 - l.Slog (^) + 1.62 {log (^) }' 



exp[-3.271396786 x 10 1 
+ 1.35365560 x 10 1 lnT c 

- 5.73932875 x 10°(lnT e ) 2 
+ 1.56315498 x 10°(lnT c ) 3 

- 2.87705600 x lO-^lnTe) 4 
+ 3.48255977 x 10- 2 (lnT c ) 5 

- 2.63197617 x 10- 3 (lnT e ) 6 
+ 1.11954395 x 10~ 4 (lnT c ) 7 
-2.03914985 x 10- 6 (lnT e ) 8 l 



1.269 x 10 



-13 (315614)1-803 



fcl2,A 

x [1.0+ (S2f2 

fc 12 , B = 2.753 xlO- 14 (^fM) 1 ' 5 °° 
x [i.o+ (iif 88 )°- 407 ]-2.242 

fc 13 = cxp[-1.801849334 x 10 1 
+ 2.36085220 x 10°lnT e 

- 2.82744300 x 10 _1 (lnT e ) 2 
+ 1.62331664 x 10- 2 (lnT c ) 3 

- 3.36501203 x 10- 2 (lnT c ) 4 
+ 1.17832978 x 10- 2 (lnT c ) 5 

- 1.65619470 x 10- 3 (lnT c ) 6 
+ 1.06827520 x 10- 4 (lnT c ) 7 

- 2.63128581 x 10- 6 (lnT e ) 8 ] 



Case A 



Case B 



9 
10 

10 

11 

12 
12 
13 



14 



14 



13 
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Table Bl continued 



14 H-+H^H + H + e" 



15 H-+H+-^H++e- 

16 He + e- -> He+ + e" 



17 He+ 



18 
19 



He + 7 



He 4 
He- 



+ H 

H+ 



20 C+ + e- 



21 0++< 



• He + H+ 

• He+ + H 

-C + 7 



+ 7 



22 C + c- ->■ C+ +c- +e" 

23 O + e- -> 0+ +e- + e" 

24 0+ + H — > O + H+ 

25 O + H+ — ► 0+ + H 

26 O + Hc+ -* 0+ + He 

27 C + H+ — > C+ + H 

28 C+ + H — > C + H+ 

29 C + He+ -> C+ + He 

30 H2 + He — + H + H + He 



31 OH + H^O + H + H 

32 HOC+ + H 2 ->■ HCO+ + H 2 

33 HOC+ + CO -» HCO+ + CO 

34 C + H2 — > CH + H 

35 CH + H -> C + H 2 



-1.78186 



fci4 = 2.5634 x lO-Te 1 

= exp[-2.0372609 x 10 1 
+ 1.13944933 x 10° lnT e 

- 1.4210135 x 10 _1 (lnT e ) 2 
+ 8.4644554 x 10- 3 (lnT e ) 3 

- 1.4327641 x 10- 3 (lnT e ) 4 
+ 2.0122503 x 10" 4 (lnT c ) 5 
+ 8.6639632 x 10- 5 (lnT c ) 6 

- 2.5850097 x 10- 5 (lnT e ) 7 
+ 2.4555012 x 10- 6 (lnT c ) 8 
-8.0683825 x 10- 8 (lnT e ) 9 ] 

fcis = 6.9 x io- 9 T-°- 35 
= 9.6 x i - 7 T-°- 90 
fc 16 = exp[-4.409864886 x 10 1 
+ 2.391596563 x 10 1 lnT c 

- 1.07532302 x lO^lnTg) 2 
+ 3.05803875 x 10°(lnT o ) 3 

- 5.6851189 x lO-^lnTe) 4 
+ 6.79539123 x 10- 2 (lnT o ) 5 
-5.0090561 x 10- 3 (lnT c ) 6 
+ 2.06723616 x 10" 4 (lnT o ) 7 
-3.64916141 x 10- 6 (lnT e ) 8 ] 

fc 17iI . r A = io-HT" 5 [12.72 - 1.615 log T 

- 0.3162(log T) 2 + 0.0493(log T) 3 ] 
fcl7,rr,B = lO" 11 ^ - 5 [11.19 - 1.676 log T 

- 0.2852(logT) 2 + 0.04433(logT) 3 ] 
fci7, di = 1.9 x 10- 3 T- 15 cxp (-41^21) 

x [i. + 0.3exp(-2^)] 
fcis = 1.25 x 10" 15 (jgj) 25 
fci 9 = 1.26 x 10- 9 T" " 



T e < 0.1 cV 



13 



' cxp 



&21 



&22 

&25 
k26 



: 4.0 x 10" 37 T 4 -' 



: 9.62 x 10" 
: 1.30 x 10" 
: 1.41 x 10" 
X exp (- 



^ 127500 ^ 



/ t 1 


-0.6 


V 300 I 


1 


( T 1 


2.49 


V 300 ) 


1 


T\ 


-1.37 



exp ^ 
exp | 



T J 
-115786.2 



10^-0.64 
10^p-0.66 



) 



+ 7.4 x io- 4 T- 15 
) [1.0 + 0.062 x exp (-i^SO 
: 6.85 x 10" 8 (0.193 + u)- 1 u°- 25 e- u 
: 3.59 x 10- 8 (0.073 + u)- 1 u°- 34 e-" 
: 4.99 x io- n T - 405 + 7.54 x lO" 10 ! 1 " - 458 
: [1.08 x io- n T°- 517 
+ 4.00 x io-K'T - 00669 ] cxp (-2I 7 ) 
= 4.991 x 10- 15 ( T ^) - 3 ^ 9 %xp (- TT ^ M ) 
-2.780 x 10" 15 (rZ^) "'^ CX P ( glioo ) 



)] 



C_i_ 

V ioooo 

16ji0.213 



fe 2 7 = 3.9 x 10 
fc 2S = 6.08 x 10~ 14 ( 
k 29 = 8.58 x io- 17 T - 757 



xl.96 / 
10000 ) CX P I 



170000 
T , 



= 3.25 x 10" 



7y0.i 



= 2.77 x lO-iST 1 - 897 
fc 30 ,i = dex [ - 27.029 + 3.801 log (T) - 29487/T] 
%)!h = dex [-2.729 - 1.75 log (T) - 23474/T] 
n CTiH e = dex [5.0792(1.0 - 1.23 x 10~ 5 (T - 2000) 
fc 3 i=6.0xl0- 9 cxp(-^0) 
fc 32 = 3.8 x 10" 10 
k 33 = 4.0 x 10" 10 
k 34 = 6.64 x 10" 
fc 35 = 1.31 x 10" 



exp 



10 exp (_80) 



T > 0.1 eV 
T < 8000 K 
T > 8000 K 



Case A 
Case B 



T < 10000 K 

T > 10000 K 

T < 7950 K 

7950 K < T < 21140 K 

T > 21140 K 

T ^ 400 K 

T > 400 K 
u = 11.26/T 



U : 



13.6/Te 



T ^ 200 K 

200 < T < 2000 K 

T > 2000 K 



15 



13 



16 
16 

17 
18 
19 

20 

21 

22 
22 
23 
24 

25 

24 
24 
26 

27 

27 
28 
29 
30 
31 
32 
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36 


CH + H 2 -> CH 2 + H 


^36 


37 


CH + C ->■ C 2 + H 


^37 


38 


CH + O -+ CO + H 


^38 


39 


CH 2 + H -> CH + H 2 


^39 


40 


CHo -4- O — > CO -4- H 4- H 


/C4Q 


41 


CH 2 + O -» CO + H 2 


&41 


42 


c 2 + o->co + c 


^42 


43 


O + H 2 OH + H 


^43 


44 


OH + H -> O + H 2 


/C44 


45 


OH + H 2 -»H 2 + H 


A,45 


46 


OH + C ->■ CO + H 




47 


OH + O -> 2 + H 


&47 


48 


OH + OH^H 2 + H 


/C48 




TT~f) 4- TT v TT^, 4- OH 


K49 


50 


i TT AIT i / — \ 

2 + H — > OH + O 


&50 


51 


2 +H 2 -> OH + OH 




52 


o 2 + c^co + o 


^52 


53 


/ ' / 1 TT , f\ T 1 T T 

OU + 11 — ► O + On 


7,, 

K53 


54 


Hj + H 2 -> H+ + H 


£,54 


55 


H+ + H -> H+ + H 2 


&55 


56 


C + H+ ->■ CH+ + H 


^56 


57 


C + H+ ->■ CH+ + H 2 


^57 


58 


C+ + H 2 -> CH+ + H 


&58 


59 


CH+ + H — > C+ + H 2 


^59 


60 


CH+ + H 2 -> CH+ + H 


^60 


61 


CH+ + O -» CO+ + H 


^61 


62 


CH 2 + H+ ->■ CH+ + H 2 


k<r,2 


63 


CH+ + H — > CH+ + H 2 


&63 


64 


CH+ + H 2 -> CH+ + H 


^64 


65 


CH+ + O -> HCO+ + H 


^65 


66 


OH+ 4- H — > CH+ 4- Ht 

V-'llg -11 > VII2 1-1 2 


^66 


67 


/~<Tt1~ 1 i" T TTPA^ 1 TT 

CH3 + O — > HCO T + H 2 


^67 


68 


o 2 + o T — » to T + 


&68 


69 


/~\+ I T T / ITT ■ 1 tt 

O t + H 2 — » OH T + H 




1 


f 1 1 Tjl , / I T T 1 1 TT 

O + rl 2 — ► Orl^ + hi 


J, 

«70 


71 
I 1 


O + rig — > Oil + ri 2 


/C7I 


70 

1 Z 


f^iTT 1 TT + * TJ~(~i+ 1 TT~ 

Orl + rig — > 11 2 <J + H 2 


k72 


71 
l o 


ott j_ . r^o+ 4_ 

Oil O — ► K^^J -f- H 


^73 


rl 


T 1 T T ' 1 TJ . TT /~i"T~ 1 TT 

Orl + H 2 — > H 2 + rl 


«74 


75 


H 2 0+ + H 2 -^H 3 0+ + H 


&76 


76 


H 2 + H+ -+H 3 0+ + H 2 


^76 


77 


H2O + C+ ->■ HCO+ +H 


fc77 


78 


H 2 + C+ -^HOC++H 


^78 


79 


H3O+ + C -> HCO+ + H 2 


^79 


80 


o 2 + C+ -r CO+ + 


^80 


81 


o 2 + C+ -» CO + O+ 


fesi 


82 


2 + CH+ -r HCO+ + OH 


^82 


83 


0+ + C -r co+ + 


/«83 


84 


CO + H+ -> HOC+ + H 2 


Al84 


85 


CO + H+ ^HCO++H 2 


&85 


86 


HCO+ + C ->■ CO + CH+ 


/«86 


87 


HCO+ + H 2 -» CO + H3O+ 


^87 



: 5.46 x 10 
: 6.59 x 10 
: 6.6 x 10" 
: 1.02 x 10 
: 6.64 x 10 
: 1.33 x 10 
: 8.0 x 10- 
: 5.0 x 10- 
5.0 x 10- 
: 3.14 x 10 
: 6.99 x 10 
: 2.05 x 10 

: 1.0 X 10- 

: 3.50 x 10 
: 1.77 x 10 

: 1.65 x 10 
: 1.59 x 10 
2.61 x 10 
: 3.16 x 10 
: 4.7 x 10- 
: 2.48 x 10 

: 1.1 X 10" 

: 2.24 x 10 
: 7.7 x 10" 
: 2.4 x 10- 
: 2.0 x 10- 

: 1.0 X 10- 

: 7.5 x 10- 

: 1.2 X 10- 

: 3.5 x 10- 

: 1.4 X 10" 
: 1.0 X 10" 
: 1.6 X 10" 

: 7.5 X 10" 

: 7.0 x 10- 

: 4.0 X 10" 

: 4.8 X 10" 

: 1.7 X 10" 
: 1.5 X 10" 

: 8.4 X 10" 

: 1.3 X 10" 

: 7.7 X 10" 

: 1.01 X 10 

: 6.4 X 10" 

: 5.9 x 10- 

: 9.0 X 10" 
: 1.8 X 10" 
: 1.0 X 10" 

: 3.8 x 10- 

: 6.2 X 10" 

: 9.1 X 10- 

: 5.2 X 10" 
: 2.7 X 10" 

: 1.7 X 10" 
: 1.1 X 10" 

: 2.5 X 10" 



-10 

-11 
11 



cxp (-±M3) 

) 



-10 cxp (_ 914 

-11 

-10 

11 



11 ( T \ C 

V 300 ) 
11 ( _T_\ 

\ 300 I 

V 300 ) _ 

(300) ex pl 




T < 2000 K 
T > 2000 K 



T ^ 300 K 
T > 300 K 



T ^ 261 K 
T > 261 K 



T ^ 295 K 
T > 295 K 



4640 \ 
T J 



cxp 



10 
9 
9 
9 
10 

i» cxp (_io^o; 

10 

10 

9 

9 

10 



33 
34 
35 
36 
37 
38 
39 
40 
41 
42 
43 
44 
34 
45 
33 
34 
46 
33 
47 
34 
33 
28 
48 
49 
28 
28 
50 
51 
51 
52 
28 
28 
53 
28 
28 
54 
28 
55 
28 
56 
28 
28 
57 
58 
59 
60 
60 
28 
53 
53 
53 
28 
61 
61 
28 
62 
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88 H 2 +He+^He + H+ fc 88 = 7.2 x 10" 15 63 

89 H 2 +Hc+ -» Hc + H + H+ fc 89 = 3.7 x lO" 14 cxp (^) 63 

90 CH + H+^CH++H fc 90 = 1.9 x 10" 9 28 

91 CH 2 + H+ -» CH+ + H fc 9 i = 1.4 x 10" 9 28 

92 CH 2 + He+ -+ C+ + He + H 2 k 92 = 7.5 x 10" 10 28 
r>_ _i_ u„+ _x n+ _i_ r 1 _i_ tio ! i « ^ 1 n— 9 28 

28 
28 
64 
65 



93 C 2 + He+ -^C+ + C + He fc 93 = 1.6 x 10" 9 

94 OH + H+^OH++H k 94 = 2.1 x 1CT 9 

95 OH + Hc+ ->■ 0+ + He + H fc 95 = 1.1 x 10" 9 

96 H 2 + H+ -> H 2 0+ + H k 96 = 6.9 x 10" 9 

97 H 2 + He+ ^ OH + He + H+ k 97 = 2.04 x 1CT 10 

98 H 2 + Hc+ -> OH+ +He + H k 9 $ = 2.86 x 10 -10 65 

65 
64 
66 

102 2 +He+ -> 0+ + + He fc 102 = 1.1 x 10" 9 66 

28 
67 

105 CO + Hc+ -^C + 0+ + Hc fei 05 = 1.4 X 10" 16 (^)"°' 5 67 

106 CO+ + H^CO + H+ fcioe = 7.5 x 10" 10 ' 68 

28 



99 H 2 + Hc+ -> H 2 0+ + He k g9 = 6.05 x 10" 11 

100 2 +H+->0+ + H fcioo = 2.0 x 10- 9 

101 2 + Hc+ -> + He fcioi = 3.3 x 10" 11 

102 2 +He+ -> 0+ + + He fc 102 = 1.1 x 10" 9 

103 0+ + C -> 2 + C+ feio3 = 5.2 x 10- 11 

104 CO + Hc+ ^C+ + + He fci04 = 1.4 x 10" 9 (^) ~°" 
"- CO + Hc+ ^C + 0+ + He k 105 = 1.4 X 10" 16 (^) "° 

CO+ + H^CO + H+ fcioe = 7.5 x 10- 10 

107 CT+H+^C + H feior = 2.3 x 10- 7 (^q)^ ' 5 

108 0~+H+-xO + H fcios = 2.3 x 10" 7 (3^)~°' 5 28 

109 Hc++H-^Hc + H fc 109 = 2.32 x 10 ~ 7 (^)~°' 5 %xp(— ^) 69 



110 H++e-^H 2 + H feno = 2.34 x lO" 8 °' 52 70 

111 H++c-^H + H + H fem = 4.36 x 10~ 8 



0.52 
-0.5 



70 
71 
72 
72 

L15 CH++e-^C + H 2 fcns =7.68 x 10- 8 (^) _0 ' 6 72 

110 CH++c-^CH 2 + H fc n6 = 7.75 x 10- 8 (t^) -0 ' 5 73 

73 
28 
74 
75 



112 CH++c-^C + H fcn2 = 7.0 x 10~ 8 (3^)"° 

113 CH++c-^CH + H fciig = 1.6 x 10- 7 (3^)"° 

114 CH++e-^C + H + H kiu = 4.03 x 10~ 7 (^)~ 
^++c-^C + H 2 feus = 7.68 x 10" 8 
3++c-^CH 2 + H fc n6 = 7.75 x 10- 8 (3^^ 

117 CH++c-^CH + H 2 fcn7 = 1.95 x 10- 7 (3^j) _0 ' 5 

118 CH+ +e- — CH + H + H fc n8 =2.0 x (^) _M 

119 OH++c-^0 + H fciig = 6.3 x 10- 9 (g^)" ' 48 

120 H 2 0+ +e" -> O + H + H k 120 = 3.05 x 10" 7 ( ^) ~°' 5 

121 H 2 0++c-^0 + H 2 fci 21 = 3.9 x 10- 8 (g^)" ' 5 75 



v 300 / 

122 H 2 0++c-^OH + H fci 22 = 8.6 x 10- 8 (g^)" ' 5 75 

123 H3O+ + e- -> H + H 2 fci 23 = 1.08 x 10~ 7 (^) ~°' 6 76 

X T X^4- . _ I V X X XX , „ „„ , „_S I T \ -0.5 



124 H3O+ + e" -x OH + H 2 k 124 = 6.02 x 10~ 8 ( ^) ~ ' 76 

125 H3O+ +e- -x OH + H + H fci 25 = 2.58 x 10~ 7 ( ^j) ~°' 6 76 

76 



126 H3O+ +c- -x + H + H 2 fci 26 =5.6 x 10- 9 (^j" ' 5 

127 Oj + c-^O + O fci 27 = 1.95 x 10~ 7 (^o) -0 ' 7 

128 CO++c--^C + fci 28 = 2.75 x 10" 7 (^)" 0M 

129 HCO+ + c" -x CO + H fci 29 = 2.76 x 10- 7 (^j) -0 ' 64 

130 HCO++c-^OH + C feiso = 2.4 x 10- 8 (^q) -0 ' 64 

131 HOC+ + e-^CO + H fei 3 i = 1.1 x 10~ 7 ( ^) _1 " 

132 r+C^CH + e" fc i32 = 1.0 x 10~ 9 

133 H-+0-+OH + e- fc 133 = 1.0 x lO" 9 28 

134 H-+OH^H 2 + e" fci 3 4 = 1.0 x 10" 10 28 

135 CT+H-xCH + c" fci 35 = 5.0 x lO" 10 

136 C-+H 2 ^CH 2 +e- k ls(i = 1.0 x lO" 13 

137 Cr+O-^CO + c" fci 37 = 5.0 x 10- 10 

138 O^+H^OH + e- fei 38 = 5.0 x lO" 10 

139 0-+H 2 ^H 2 + e" fc 139 = T.O x lO" 10 28 

140 0~+C->CO + c- /C140 = 5.0 x 10- 10 28 

141 H 2 +H+^H++7 fci4i = 1.0 x 10- 16 



77 

78 

79 

79 

28 
28 



28 
28 
28 
28 
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Table Bl continued 



142 C + e-^CT+7 fci42 = 2.25 X 10~ 15 81 

143 C + H^CH + 7 fci43 = 1.0 x 1CT 17 82 

144 C + H 2 ^CH 2 +7 fci44 = 1.0 x 10~ 17 82 

145 C + C^C 2 +7 fews = 4.36 x 10" 18 (3^)°' 35 exp (-i^) 83 

146 C + O^CO + 7 fci4 6 = 2.1 x 10~ 19 T ^ 300 K 84 

= 3.09xl0- 17 (^)°' 33 exp(-±f^) T> 300K 85 

147 C+ + H-^CH++7 fci47 = 4.46 x 10- 16 T-°' 5 exp (- 86 



-0.2 



148 C+ + H 2 ^CH+ + 7 fc 148 = 4.0 x 10" 16 (^q) ' 87 

149 C+ + O^CO++7 fci4 9 = 2.5 x 10" 18 T ^ 300 K 84 

= 3.14 x 10~ 18 (g^j)" 15 exp (^) T> 300K 

150 + e~^0-+7 fci 50 = 1.5 x 10~ 15 28 

151 O + H^OH + 7 fcisi = 9.9 x 10- 19 (g^j) -0 ' 38 28 

152 + 0^0 2 +7 fci 52 = 4.9 x 10" 20 (g^) 1 ' 58 82 

153 OH + H^H 2 + 7 k 153 = 5.26 x 10" 18 ( ^) ~ 5 ' 22 exp (-f ) 88 

154 H + H + H^H 2 +H fci 54 = 1.32 x 10~ 32 ( ^) ~°' 38 T ^ 300 K 89 

= 1.32 x 10" 32 (g^o) -10 T> 300K 90 

155 H + H + H 2 ~»H 2 + H 2 k 155 = 2.8 x io- 31 T-°- 6 91 

156 H + H + He^H 2 + Hc fc 156 = 6.9 x icr^T" - 4 92 

157 C + C + M^C 2 +M fci57 = 5.99 x 10~ 33 (g^)" 1 ' 6 T ^ 5000 K 93 

= 5.99 x 10- 33 (g^) -0 ' 64 exp (^5) T> 5000K 94 

158 C + O + M^CO + M fci 58 = 6.16 x 10" 29 (g^)" 3 ' 08 T ^ 2000 K 35 

= 2.14 x 10- 29 (ggj)" 3 ' 08 exp (3«4) T> 2000K 67 

159 C+ + O + M — > CO+ + M fci 59 = 100 x k 210 67 

160 C + 0+ + M -> CO+ + M fcieo = 100 x k 2 io 67 

161 O + H + M^OH + M fci 6 i = 4.33 x 10" 32 (g^)" 10 43 

162 OH + H + M->H 2 + M fc 162 = 2.56 x 10" 31 ( ^) ~ 2 '° 35 

163 + + M^0 2 + M fci 63 = 9.2 x 10~ 34 (g^j) -10 37 

164 O + CH -» HCO+ + e- fci 64 = 2.0 x 10~ n (g^) ' 44 95 

165 H + H(s)^H 2 fci 65 = 3.0 x 10- 18 T°- 5 /A[1.0 + 0.04(T + T d )°- 5 / A = [l.O + 10 4 exp (-f^)] _1 96 

+ 0.002T + 8 x lO^T 2 ]- 1 



Notes: T and T c are the gas temperature in units of Kelvin and eV respectively, while is the dust temperature in Kelvin 
H(s) denotes a h ydrogen a tom a dsor bed on a grain surfac e. References are to the prim ary source of data for each reacti on 
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Table B2. List of photochemical reactions included in our chemical model 



No. Reaction Optically thin rate (s 1 ) 7 Ref. 



166 


H- 4 


-7^H + e" 


K166 = 7.1 


X 


10" 


■7 


0.5 


1 


167 


H+ 4 


-7^H + H+ 


Rl67 = 1.1 


X 


10" 


-9 


1.9 


2 


168 


H 2 + 


7^H + H 


R168 = 5.6 


X 


10" 


■'11 


See 


3 


169 


H+4 


7 -+ H 2 + H+ 


Rl69 = 4.9 


X 


10- 


-13 


1.8 


4 


170 


H+4 


7 -> H+ 4- H 


Rno = 4.9 


X 


10- 


■13 


2.3 


4 


171 


C 4-7^0+4- e" 


Hl71 =3.1 


X 


10- 


■10 


3.0 


5 


172 


C" 4 


-7 -» C + e" 


R172 = 2.4 


X 


10" 


■7 


0.9 


6 


173 


CH 4 


-7^C4H 


i?173 = 8.7 


X 


10" 


-10 


1.2 


7 


174 


CH 4 


- 7 -> CH+ +e" 


R 174 = 7.7 


X 


10- 


-10 


2.8 


8 


175 


CH+ 


+ 7^C + H+ 


Rm = 2.6 


X 


10- 


■10 


2.5 


7 


176 


CH 2 


4- 7 -> CH + H 


i?176 = 7.1 


X 


10- 


■10 


1.7 


7 


177 


CH 2 


4- 7 -+ CH+ 4- e" 


iim = 5.9 


X 


10- 


■10 


2.3 


6 


178 


CH+ 


+ 7 -» CH+ 4- H 


#178 = 4.6 


X 


10- 


■10 


1.7 


9 


179 


CH+ 


+ 7 -> CH+ 4- H 


#179 = 1.0 


X 


10- 


9 


1.7 


6 


180 


CH 3 + 


+ 7 -» CH+ 4- H 2 


#180 = 1.0 


X 


10- 


9 


1.7 


6 


181 


C 2 4- 


7 -> C + C 


#181 = 1.5 


X 


10- 


■10 


2.1 


7 


182 


O- 4 


-7 -> O + e~ 


#182 = 2.4 


X 


10- 


■7 


0.5 


6 


18.3 


OH 4 


-7^0 + H 


#183 = 3.7 


X 


10- 


■10 


1.7 


10 


184 


OH 4 


-7-» OH+ +e" 


#184 = 1.6 


X 


10- 


■12 


3.1 


6 


1 C K 

loO 


OH+ 


+ 7 -> 4-H+ 


d in 

-K185 — t.u 


X 


1U 


■1 2 


1.0 


1 


186 


H 2 


+ 7 -» OH + H 


#186 = 6.0 


X 


10- 


■10 


1.7 


11 


187 


H 2 


+ 7 -> H 2 0+ +e~ 


#187 = 3.2 


X 


10- 


■11 


3.9 


8 


188 


H 2 0" 


h 4-7^H+4-0 


#188 = 5.0 


X 


10- 


■11 


See 3T2] 


12 


189 


H 2 0" 


h 4-7 -» H+ 4- OH 


#189 = 5.0 


X 


10- 


■11 


See 322] 


12 


190 


H 2 0" 


h 4-7 -» 0+ 4-H 2 


#190 = 5.0 


X 


10- 


-11 


See Eg] 


12 


191 


H 2 0" 


h 4-7 -> OH+ +H 


#191 = 1.5 


X 


10- 


■10 


See 322] 


12 


192 


H3O- 


h 4-7 -> H+ 4-H 2 


#i 92 = 2.5 


X 


10- 


■11 


See |22] 


12 


19.3 


H3O- 


h 4- 7^ H+4- OH 


#193 = 2.5 


X 


10- 


■11 


See EH 


12 


194 


H3O- 


h 4-7^H 2 0+ + H 


#194 = 7.5 


X 


10- 


■1 2 


See EH 


12 


195 


H3O- 


h 4- 7 ~» OH+ + H 2 


#195 = 2.5 


X 


10- 


■11 


See 3221 


12 


196 


O2 + 


7 -» 0+ 4-e" 


#196 = 5.6 


X 


10- 


■11 


3.7 


7 


197 


2 + 


7^0 + 


#197 = 7.0 


X 


10- 


■10 


1.8 


7 


198 


CO 4 


-7-*C4-0 


#198 = 2.0 


X 


10- 


-10 


See 3221 


13 



Note: Rates are computed assuming the standard interstellar radiation field from lDrainel lll978l) . with field strength Go = 1.7 in lHabinel 
units. 7 quantifies the dependence of the rate on the visual extinction Ay in optically thick gas: # t hick = #;hin ex P( — 7^4-v) ( see 

Eq.Q. 

Refer ence s: 1: Ide Jond dl972h. 2: IPunnl dl968t). 3: iDraine fc Bertoldil lll99rj). 4: Ivan Dishoeckl Jl987h. 5: IVerner et al.l 
dl996h. 6: iLe Teuff. Millar fc MarkwicM jl^OoT) . 7: iRoberee et al.l lll99lh. 8: Ivan Dishoeckl dl988h. 9: Ivan DishoeclJ 1 I2OO6I) . 10: 
Ivan Dishoeck fc Dalgarnol l|l984h . ll:lLed l|l984l ), 12: Sternberg &c Dalgarn ol <1995h . 13: Ivan Dishoeck fc Blackl l|l988l ) 
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Table B3. List of reactions included in our chemical model that involve cosmic rays or cosmic-ray induced UV emission 



No. 


Reaction 




Rate (s -1 ^ 1 ) 


Ref. 


199 


H + c.r. — » 


H+ + c~ 


^199 


= 1.0 




200 


He -h c.r. — 


♦ He+ + e~ 


R200 


= 1.1 


1 


201 


H 2 + c.r. - 


+ H+ + H + c~ 


R201 


= 0.037 


1 


202 


H2 + c.r. — 


* H + H 


R202 


= 0.22 


1 


203 


H 2 + c.r. - 


+ H+ +H _ 


R203 


= 6.5 x 10~ 4 


1 


204 


H 2 + c.r. — 


» h+ 4. e - 


R204 


= 2.0 


1 


205 


C + c.r. -> 


C+ +c- 


R205 


— 0.0 


1 


206 


O + c.r. -* 


0+ +c- 


R2O6 


= 5.7 


1 


207 


CO + c.r. - 


-> CO+ + e" 


R207 


= 6.5 


1 


208 


C + 7c. r . - 


> C+ +e- 


R20S, 


= 2800 


2 


209 


CH + 7 c.r. 


^C + H 


R209 


= 4000 


3 


210 


CH+ + 7c ., 


■. -» C+ + H 


R210 


= 960 


3 


211 


CH 2 + 7c.r 


. -» CH+ + e- 


R211 


= 2700 


1 


212 


CH 2 + 7c.r 


. -> CH + H 


R212 


= 2700 


1 


213 


C 2 + 7c. r. - 


^C + C 


R213 


= 1300 


3 


214 


OH + 7c.r. 


^O + H 


R214 


= 2800 


3 


215 


H 2 + 7 c . r 


. -» OH + H 


R215 


= 5300 


3 


216 


2 + 7c. r. - 


-^O + O 


R216 


= 4100 


3 


217 


2 + 7c. r. - 


-O+4-e" 


R217 


= 640 


3 


218 


CO + 7c.r. 


^C + O 


R218 


= 0.21T 1 / 2 ccH 2 a;co' 2 


1 



Note: Rates are quoted relative to the cosmic ray ionization rate of atomic hydrogen, which is an adjustable parameter in our 
models. Rates for cosmic-ray induced photoionizations and photodissociations (reactions 208-218) arc quoted assuming a grain albedo 
lu = 0.6, following reference 1. In the expression for -R 2 i8, T is the gas temperature in Kelvin and xn 2 and xqq are the fractional 

abundances of H 2 a nd CO, respectively. 

References: 1: ILe Teuff. Millar fc Markwickl (|200(Jl . 2: iGredel. Lepp fc Dalgarnol l|l987l l. 3: iGredel et all [|l989l l. 4: 
iMalonev. Ho llcnbach &: Tielen j dl996h 



